Package stallone.datasequence.io

Source Code of stallone.datasequence.io.PairIterable

package stallone.datasequence.io;

import static stallone.api.API.*;

import java.io.IOException;
import java.nio.ByteBuffer;
import java.nio.ByteOrder;
import java.nio.FloatBuffer;
import java.util.BitSet;
import java.util.Iterator;
import java.util.logging.Level;
import java.util.logging.Logger;

import stallone.api.datasequence.DataSequence;
import stallone.api.datasequence.IDataList;
import stallone.api.datasequence.IDataReader;
import stallone.api.datasequence.IDataSequence;
import stallone.api.doubles.Doubles;
import stallone.api.doubles.IDoubleArray;
import stallone.io.NicelyCachedRandomAccessFile;

/**
* Reader for charmm / xplor / namd dcd files.
*
* <p>File specification and ideas for source code gained from sources of MolTools package (Frank Noe). source of VDM
* plugin: dcdplugin.c ( http://www.ks.uiuc.edu/Research/vmd/plugins/doxygen/dcdplugin_8c-source.html )</p>
*
* @author  Martin Senne
*/
public class DcdReader implements IDataReader
{
    private String filename;
    private boolean format__dcd_is_charmm;
    private boolean format__dcd_has_4dims;
    private boolean format__dcd_has_extra_block;
    private ByteOrder byteOrder;
    private NicelyCachedRandomAccessFile niceRandomAccessFile;
    /** Number of sets of coordinates (NSET), that is number of frames. */
    private int numberOfFrames;
    /** Starting timestep (ISTART). */
    private int startingTimestep;
    /** Store NSAVC, the number of timesteps between dcd saves (NSAVC. */
    private int numberOfTimestepsBetweenSaves;
    /** Number of atoms. */
    private int numberOfAtoms;
    /** Number of fixed atoms (NAMNF). */
    private int numberOfFixedAtoms;
    /** Timestep delta (DELTA). */
    private double delta;
    private long frameAreaStartingPosition;
    private long framesizeRegular;
    private long framesizeWithFixed;
    private int nextFrameIndex;

    /** Selection for reading */
    private int[] selection;
    private BitSet selected;
   
    private Logger logger;

    /**
     * empty DcdReader, use {@link DcdReader#setSource(String)},
     * {@link DcdReader#open()} and {@link DcdReader#scan()}
     */
    public DcdReader()
    {
    }
   
    /**
     * constructs DcdReader with given filename
     * @param _filename
     * @throws IOException
     */
    public DcdReader(String _filename) throws IOException
    {
        this.filename = _filename;
        this.logger = Logger.getLogger(DcdReader.class.getName());
        initialize();
    }

    @Override
    public void setSource(String name)
    {
        this.filename = name;
    }
   
    @Override
    public void scan()
            throws IOException
    {
        this.initialize();
    }

    private void initialize() throws IOException
    {
        niceRandomAccessFile = new NicelyCachedRandomAccessFile(this.filename);
        byteOrder = detectEndianess();
        readHeader();
        niceRandomAccessFile.changePageSize((int) framesizeRegular);
       
        selection = intArrays.range(numberOfAtoms);
        selected = new BitSet(numberOfAtoms);
        selected.set(0, numberOfAtoms);
    }

    /**
     * First record length indicator (84 bytes) gives the clue of which endian byte order is used 54 00 00 00 =>
     * ByteOrder.LITTLE_ENDIAN 00 00 00 54 => ByteOrder.BIG_ENDIAN.
     *
     * @return  appropiate byte order
     */
    private ByteOrder detectEndianess() throws IOException
    {
        niceRandomAccessFile.seek(0);

        ByteBuffer recordLengthBuffer = niceRandomAccessFile.readToBuffer(4);

        niceRandomAccessFile.seek(0);

        if (recordLengthBuffer.order(ByteOrder.LITTLE_ENDIAN).getInt() == 84)
        {
            logger.log(Level.FINE, "Little endian encoding.");
            return ByteOrder.LITTLE_ENDIAN;
        }
        else if (recordLengthBuffer.order(ByteOrder.BIG_ENDIAN).getInt() == 84)
        {
            logger.log(Level.FINE, "Big endian encoding.");
            return ByteOrder.BIG_ENDIAN;
        }
        else
        {
            throw new IOException("Invalid record length of first record.");
        }
    }

    /**
     * Read next record.
     *
     * @param   buffer  is the target
     *
     * @return
     *
     * @throws  IOException
     */
    private ByteBuffer readNextRecord() throws IOException
    {
        // read record length (int)
        int recordLength = niceRandomAccessFile.readToBuffer(4).order(byteOrder).getInt();

        // make sure, data fits in one buffer
        // niceRandomAccessFile.fitsInBuffer( recordLength );

        ByteBuffer recordByteBuffer = niceRandomAccessFile.readToBuffer(recordLength).order(byteOrder);

        // read away second length indicator at the end of the record
        int verifyRecordLength = niceRandomAccessFile.readToBuffer(4).order(byteOrder).getInt();

        // check record length
        if (recordLength == verifyRecordLength)
        {
            logger.log(Level.FINEST, "Record (length=" + recordLength + ") found and read.");
            return recordByteBuffer;
        }
        else
        {
            throw new IOException("Invalid file data.");
        }
    }

    /**
     * Skip record.
     *
     * @throws  IOException
     */
    private void skipNextRecord() throws IOException
    {

        // read record length (int)
        int recordLength = niceRandomAccessFile.readToBuffer(4).order(byteOrder).getInt();

        niceRandomAccessFile.skipBytes(recordLength);

        int verifyRecordLength = niceRandomAccessFile.readToBuffer(4).order(byteOrder).getInt();

        // check record length
        if (recordLength == verifyRecordLength)
        {
            logger.log(Level.FINEST, "Record (length=" + recordLength + ") found and skipped.");
            return;
        }
        else
        {
            throw new IOException("Invalid file data.");
        }
    }

    /**
     * Reads the header information. This are the first 3 records: 1. Main header record 2. Description record 3. Number
     * of atoms record
     */
    private void readHeader() throws IOException
    {

        // unset all flags
        format__dcd_is_charmm = false;
        format__dcd_has_4dims = false;
        format__dcd_has_extra_block = false;

        ByteBuffer mainHeaderRecord = readNextRecord();

        // check for proper size of first record and CORD identifier
        if ((mainHeaderRecord.remaining() == 84)
                && (mainHeaderRecord.get(0) == 'C')
                && (mainHeaderRecord.get(1) == 'O')
                && (mainHeaderRecord.get(2) == 'R')
                && (mainHeaderRecord.get(3) == 'D'))
        { // this is cord header

            // CHARMm-genereate DCD files set the last integer in the
            // header, which is unused by X-PLOR, to its version number.
            // Checking if this is nonzero tells us this is a CHARMm file
            // and to look for other CHARMm flags.
            if (mainHeaderRecord.getInt(80) != 0)
            {
                format__dcd_is_charmm = true;

                logger.log(Level.FINE, "CHARMM format DCD file (also NAMD 2.1 and later).");

                if (mainHeaderRecord.getInt(44) != 0)
                {
                    format__dcd_has_extra_block = true;

                    logger.log(Level.FINE, "  Has extra block.");
                }

                if (mainHeaderRecord.getInt(48) != 0)
                {
                    format__dcd_has_4dims = true;

                    logger.log(Level.FINE, "  Has 4-dims.");
                }
            }
            else
            {
//                format__dcd_is_xplor = true; // must be an X-PLOR format DCD file
                logger.log(Level.FINE, "X-PLOR format DCD file (also NAMD 2.0 and earlier).");
            } // end if-else

            // Store the number of sets of coordinates (NSET)
            this.numberOfFrames = mainHeaderRecord.getInt(4);

            logger.log(Level.FINE, "  Number of sets of coordinates: " + this.numberOfFrames);

            // Store ISTART, the starting timestep
            this.startingTimestep = mainHeaderRecord.getInt(8);

            logger.log(Level.FINE, "  Starting timestep            : " + this.startingTimestep);

            // Store NSAVC, the number of timesteps between dcd saves
            this.numberOfTimestepsBetweenSaves = mainHeaderRecord.getInt(12);

            logger.log(Level.FINE, "  Timesteps between dcd saves  : " + this.numberOfTimestepsBetweenSaves);

            // Store NAMNF, the number of fixed atoms */
            this.numberOfFixedAtoms = mainHeaderRecord.getInt(36);

            logger.log(Level.FINE, "  Number of fixed atoms        : " + this.numberOfFixedAtoms);

            if (this.numberOfFixedAtoms > 0)
            {
                throw new RuntimeException("Sorry, read of fixed atom coordinates currently not implemented.");
            }

            // Timestep delta. Read in the timestep, DELTA
            // Note: DELTA is stored as a double with X-PLOR but as a float with CHARMm
            if (format__dcd_is_charmm)
            {
                float ftmp = mainHeaderRecord.getFloat(40); // is this safe on Alpha?
                this.delta = (double) ftmp;
            }
            else
            {
                this.delta = mainHeaderRecord.getDouble(40);
            }

            logger.log(Level.FINE, "  Timestep                     : " + this.delta);

            ByteBuffer descriptionHeaderRecord = readNextRecord();
            int lines = descriptionHeaderRecord.getInt();
            byte[] line = new byte[80];

            for (int i = 0; i < lines; i++)
            {
                descriptionHeaderRecord.get(line);
            }

            ByteBuffer noOfAtomsRecord = readNextRecord();
            this.numberOfAtoms = noOfAtomsRecord.getInt(0);

            logger.log(Level.FINE, "No of atoms : " + this.numberOfAtoms);

            // get position framestart
            this.frameAreaStartingPosition = niceRandomAccessFile.getFilePointer();

            // calculate framesize
            this.framesizeRegular = 0;

            long dims = 3;

            if (format__dcd_has_4dims)
            {
                dims = 4;
            }

            long fullExtraframesize = 0;

            if (format__dcd_has_extra_block)
            {
                fullExtraframesize = 4 + 48 + 4; // length indicator + extra block (unit cell) + length indicator
            }

            // Size of a regular frame (every frame, if no fixed atoms.
            // Or every frame expect the first when fixed coordinates are used.
            long cordRecordSize = (numberOfAtoms - numberOfFixedAtoms) * 4;
            long fullRecordSize = 4 + cordRecordSize + 4;

            // First frame size, if fixed coordinates are used
            long cordRecordSizeWithFixed = numberOfAtoms * 4;
            long fullRecordSizeWithFixed = 4 + cordRecordSizeWithFixed + 4;

            this.framesizeWithFixed = ((fullRecordSizeWithFixed * dims) + fullExtraframesize);
            this.framesizeRegular = ((fullRecordSize * dims) + fullExtraframesize);
            this.nextFrameIndex = 0;

            // integrity check by file length
            long endpos = frameAreaStartingPosition + (numberOfFrames * framesizeRegular);

            if (niceRandomAccessFile.length() == endpos)
            {
                logger.log(Level.FINE, "Checking file size vs. calculated number of frames .... OK.");
            }
            else
            {
                throw new IOException("Real file size does not match the calculated file size. "
                        + "Maybe an interrupt has occured during creation of trajectory file." + "Real: "
                        + niceRandomAccessFile.length() + "  Expected: " + endpos);
            }

        }
        else
        { // not the cord header ..... => error
            throw new IOException("Invalid file format. Header corrupted or not found.");
        } // end if-else
    }

    /**
     * Read a frame at current position (advances in file).
     *
     * @param   factory  to construct vector
     *
     * @return  vector with frame data, that is
     *
     *          <pre>x_0, y_0, z_0, x_1, y_1, z_1, ....</pre>
     *
     * @throws  IOException
     */
    private IDoubleArray readNextFrame(IDoubleArray target) throws IOException
    {
      // ensure frame fits into buffer.
        niceRandomAccessFile.fitsInBuffer((int) framesizeRegular);

        if (format__dcd_has_extra_block)
        {
            // read away extra block, which contains unit cell per frame
            skipNextRecord();
        }

        int n = this.numberOfAtoms - this.numberOfFixedAtoms;

        FloatBuffer xBuffer = readNextRecord().asFloatBuffer();
        FloatBuffer yBuffer = readNextRecord().asFloatBuffer();
        FloatBuffer zBuffer = readNextRecord().asFloatBuffer();

        int r = 0;
        float x,y,z;
        for (int i = 0; i < n; i++)
        {
            x = xBuffer.get();
            y = yBuffer.get();
            z = zBuffer.get();
            if (selected.get(i))
            {
                target.set(r, 0, x);
                target.set(r, 1, y);
                target.set(r, 2, z);
                r++;
            }
        }

        // increment last position
        nextFrameIndex++;

        if (format__dcd_has_4dims)
        {
            skipNextRecord();
        }

        return target;
    }

    /**
     * Get position in file for frame with index frameIndex.
     *
     * @param   frameIndex
     *
     * @return  position.
     */
    private long getPositionOfFrame(int frameIndex)
    {

        if (frameIndex > 0)
        {
            return frameAreaStartingPosition + framesizeWithFixed + ((frameIndex - 1) * framesizeRegular);
        }
        else if (frameIndex == 0)
        {
            return frameAreaStartingPosition;
        }
        else
        {
            throw new IndexOutOfBoundsException("Index of frame is negative.");
        }
    }

    @Override
    public int size()
    {
        return this.numberOfFrames;
    }

    @Override
    public int dimension()
    {
        return this.numberOfAtoms * 3;
    }

    public double getTimeOfFrame(int frameIndex)
    {
        int currentStep = startingTimestep + numberOfTimestepsBetweenSaves * frameIndex;
        return ((double) currentStep) * delta;
    }

    public IDoubleArray get(int frameIndex, IDoubleArray target)
    {
        try
        {
            // position if necessary
            if ((frameIndex < 0) || (this.numberOfFrames <= frameIndex))
            {
                throw new IndexOutOfBoundsException("Invalid frame index " + frameIndex + ".");
            }

            if (frameIndex != this.nextFrameIndex)
            {
                this.nextFrameIndex = frameIndex;
                niceRandomAccessFile.seek(getPositionOfFrame(frameIndex));

                logger.log(Level.FINEST, "Repositioning for frame " + frameIndex + " to "
                            + niceRandomAccessFile.getFilePointer());
            }
            else
            {
                logger.log(Level.FINEST, "File position is okay. Frame " + frameIndex + " is at "
                            + niceRandomAccessFile.getFilePointer());
            }

            return readNextFrame(target);
        }
        catch (IOException ex)
        {
            logger.log(Level.SEVERE, null, ex);
            throw new RuntimeException("Problems reading DCD, caught I/O exception."
                    + " Message: " + ex.getMessage());
        } // end try-catch
    }

    @Override
    public void close() throws IOException
    {
        niceRandomAccessFile.close();
    }

    @Override
    public void open() throws IOException
    {
        niceRandomAccessFile.open();
    }

    @Override
    public double getTime(int frameIndex)
    {
        return delta*frameIndex;
    }

    @Override
    public void select(int[] _selection)
    {
        if (_selection == null)
            _selection = intArrays.range(numberOfAtoms);
        this.selection = _selection;
       
        // clear current selection
        selected.clear();
        // adopt to given selection
        for (int i=0; i<_selection.length; i++)
            selected.set(_selection[i]);
    }
   
    @Override
    public int[] getSelection()
    {
        return this.selection;
    }

   
    @Override
    public IDoubleArray get(int frameIndex)
    {
        return get(frameIndex, Doubles.create.array(selection.length,3));
    }

    @Override
    public IDoubleArray getView(int index)
    {
        return get(index);
    }

    @Override
    public Iterator<IDoubleArray> iterator()
    {
        return new DataReaderIterator(this);
    }

    @Override
    public Iterator<IDoubleArray[]> pairIterator(int spacing)
    {
        return new DataReaderPairIterator(this, spacing);
    }

    @Override
    public Iterable<IDoubleArray[]> pairs(int spacing)
    {
        class PairIterable implements Iterable<IDoubleArray[]>
        {
            private IDataReader seq;
            private int spacing = 1;

            public PairIterable(IDataReader _seq, int _spacing)
            {
                this.seq = _seq;
                this.spacing = _spacing;
            }

            @Override
            public Iterator<IDoubleArray[]> iterator()
            {
                return (new DataReaderPairIterator(seq, spacing));
            }
        }
        return new PairIterable(this,spacing);
    }
   
    @Override
    public IDataSequence load()
    {
        IDataList res = DataSequence.create.list();
        for (int i=0; i<this.size(); i++)
            res.add(get(i));
        return res;
    }

    public String getFileName()
    {
        return filename;
    }

    @Override
    public long memorySize()
    {
        return this.size()*this.dimension()*8;
    }
}
TOP

Related Classes of stallone.datasequence.io.PairIterable

TOP
Copyright © 2018 www.massapi.com. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.