text:
function aedat = ImportAedatDataVersion3(aedat)
%{
This is a sub-function of importAedat - it process the data where the aedat
file format is determined to be 3.x
The .aedat file format is documented here:
http://inilabs.com/support/software/fileformat/
This function is based on a combination of the loadaerdat function and
sensor-specific address interpretation functions.
There is a single input "info", a structure with the following
fields:
	- beginningOfDataPointer - Points to the byte before the beginning of the
		data
	- fileHandle - handle of the open aedat file in question
	- fileFormat - indicates the version of the aedat file format - should
		be 3.xx
	- startTime (optional) - if provided, any data with a timeStamp lower
		than this will not be returned.
	- endTime (optional) - if provided, any data with a timeStamp higher than
		this time will not be returned.
	- startPacket (optional) Any packets with a lower count that this will not be returned.
	- endPacket (optional) Any packets with a higher count that this will not be returned.
	- modPacket (optional) For sparse sampling of a file,
        packets are imported if mod(packetNumber,modPacket) = 0
	- source - a string containing the name of the chip class from
		which the data came. Options are:
		- dvs128
		- davis240a
		- davis240b
		- davis240c
		- davis128mono
		- davis128rgb
		- davis208rgbw
		- davis208mono
		- davis346rgb
		- davis346mono
		- davis346bsi
		- davis640rgb
		- davis640mono
		- hdavis640mono
		- hdavis640rgbw
		- das1
		- file
		- network
	- dataTypes (optional) cellarray. If present, only data types specified
		in this cell array are returned. Options are:
		special; polarity; frame; imu6; imu9; sample; ear; config; 1dPoint;
		2dPoint; 3dPoint 4dPoint; dynapse

The output is a structure with 2 fields.
	- info - the input structure, embelished with other data
	- data a structure which contains one structure for each type of
		data present. These structures are named according to the type of
		data; the options are:
		- special
		- polarity
		- frame
		- imu6
		- sample
		- ear
		Other data types supported in aedat3.0 are not implemented
		because no chip class currently implements them.
		Within each of these structures, there are typically a set of column
			vectors containing timeStamp, a valid bit and then other data fields,
			where each vector has the same number of elements.
			There are some exceptionss to this.
			In detail the contents of these structures are:
		- special
			- valid (colvector bool)
			- timeStamp (colvector uint64)
			- address (colvector uint32)
		- polarity
			- valid (colvector bool)
			- timeStamp (colvector uint64)
			- x (colvector uint16)
			- y (colvector uint16)
			- polarity (colvector bool)
		- frame
			- valid (bool)
			- roiId (uint8)
			- colorChannel (uint8)
			- colorFilter (uint8)
			- frame timeStamp start (uint64)
			- frame timeStamp end (uint64)
			- timeStampExposureStart (uint64)
			- timeStampExposureEnd (uint64)
			- samples (cellArray, with one cell for each frame; cells
				contain a matrix of uint16 row*col*chn, where row is the number of rows,
				col is the number of columns, and chn is the number of
				(colour) channels. Where frames have only one channel, the
				third dimension is squeezed out.
			- xLength (uint32)
			- yLength (uint32)
			- xPosition (uint32)
			- yPosition (uint32)
		- imu6
			- valid (colvector bool)
			- timeStamp (colvector uint64)
			- accelX (colvector single)
			- accelY (colvector single)
			- accelZ (colvector single)
			- gyroX (colvector single)
			- gyroY (colvector single)
			- gyroZ (colvector single)
			- temperature (colvector single)
		- sample
			- valid (colvector bool)
			- timeStamp (colvector uint64)
			- sampleType (colvector uint8)
			- sample (colvector uint32)
		- ear
			- valid (colvector bool)
			- timeStamp (colvector uint64)
			- position (colvector uint8)
			- channel (colvector uint16)
			- neuron (colvector uint8)
			- filter (colvector uint8)

Implementation: There is an efficient implementation of startPacket and
EndPacket, since the correct file locations to read from can be determined
in advance.
There are two possibilities for handling startTime and endTime; one is with
strict capture of events within the prescribed time window. The other is a
loose interpretation with capture of all events whose packets start within
the prescribed time window. It is much more efficient to implement this
second approach, and nevertheless allows a calling function to iterate
through all the data in bite-sized chunks.
There is a switch info.strictTime - if this is present and
true then the strict time approach is used, otherwise the packet-based time
approach is used. In the strict approach, data must be accumulated from the
beginning of the file and then cut off once packets' first timestamp is
greater than endTime. In the strict approach, frames are considered part of
the time window if the time which is half way between exposure start and
exposure end is within the time window.
2016_06_27 Strict time handling is not implemented yet.
Since there is no way to know in advance how big the data vectors must be,
the data vectors/matrices are started off when they are needed
and are grown by a factor of 2 each time they need to be enlarged.
At the end of the run they are clipped to the correct size.
%}

dbstop if error

info = aedat.info;
importParams = aedat.importParams;
fileHandle = importParams.fileHandle;

% Check the startEvent and endEvent parameters
if isfield(importParams, 'startPacket')
    startPacket = importParams.startPacket;
else
	startPacket = 1;
end

if isfield(importParams, 'endPacket')
	endPacket = importParams.endPacket;
else
	endPacket = inf;
end

if startPacket > endPacket
	error([	'The startPacket parameter is ' num2str(startPacket) ...
		', but the endPacket parameter is ' num2str(endPacket) ]);
end
if isfield(importParams, 'startEvent')
	error('The startEvent parameter is set, but range by events is not available for .aedat version 3.x files')
end
if isfield(importParams, 'endEvent')
	error('The endEvent parameter is set, but range by events is not available for .aedat version 3.x files')
end
if isfield(importParams, 'startTime')
    startTime = importParams.startTime;
else
	startTime = 0;
end
if isfield(importParams, 'endTime')
    endTime = importParams.endTime;
else
	endTime = inf;
end
if startTime > endTime
	error([	'The startTime parameter is ' num2str(startTime) ...
		', but the endTime parameter is ' num2str(endTime) ]);
end
% By default, throw away timeStampFrameStart/End,
% renaming timeStampExposureStart/End to timeStampStart/End
if isfield(importParams, 'simplifyFrameTimeStamps')
  simplifyFrameTimeStamps = importParams.simplifyFrameTimeStamps;
else
  simplifyFrameTimeStamps = true;
end

% By default, throw away the valid flags,
% and any events which are set as invalid.
if isfield(importParams, 'validOnly')
  validOnly = importParams.validOnly;
else
  validOnly = true;
end

% By default, do not skip any packets.
if isfield(importParams, 'modPacket')
  modPacket = importParams.modPacket;
else
  modPacket = 1;
end

%By default, import the full data, rather than just indexing the packets
if isfield(importParams, 'noData')
  noData = importParams.noData;
else
  noData = false;
end

% By default, import all data types
if isfield(importParams, 'dataTypes')
    allDataTypes = false;
    dataTypes = importParams.dataTypes;
else
    allDataTypes = true;
end

packetCount = 0;

% Has this file already been indexed in a previous pass?
if isfield(info, 'packetPointers')
	packetTypes = info.packetTypes;
	packetPointers = info.packetPointers;
	packetTimeStamps = info.packetTimeStamps;
elseif endPacket < inf
    packetTypes = ones(endPacket, 1, 'uint16');
    packetPointers = zeros(endPacket, 1, 'uint64');
    packetTimeStamps = zeros(endPacket, 1, 'uint64');
else
    packetTypes = ones(1000, 1, 'uint16');
    packetPointers = zeros(1000, 1, 'uint64');
    packetTimeStamps = zeros(1000, 1, 'uint64');
end

if noData == false
    % Create structures to hold the output data

    specialNumEvents	= 0;
    specialValid		= false(0); % initialising this tells the first pass
                                     % to set up the arrays with the size
                                     % necessary for the initial packet
    specialDataMask = hex2dec('7E');
    specialDataShiftBits = 1;

    polarityNumEvents	= 0;
    polarityValid		= false(0);
    polarityYMask = hex2dec('1FFFC');
    polarityYShiftBits = 2;
    polarityXMask = hex2dec('FFFE0000');
    polarityXShiftBits = 17;

    frameNumEvents	= 0;
    frameValid		= false(0);
    frameColorChannelsMask = hex2dec('E');
    frameColorChannelsShiftBits = 1;
    frameColorFilterMask = hex2dec('70');
    frameColorFilterShiftBits = 4;
    frameRoiIdMask = hex2dec('3F80');
    frameRoiIdShiftBits = 7;

    imu6NumEvents	= 0;
    imu6Valid		= false(0);

    sampleNumEvents	= 0;
    sampleValid		= false(0);

    earNumEvents	= 0;
    earValid		= false(0);

    point1DNumEvents = 0;
    point1DValid	= false(0);
    point1DTypeMask = hex2dec('FE');
    point1DTypeShiftBits = 1;

    point2DNumEvents = 0;
    point2DValid	= false(0);
    point2DTypeMask = hex2dec('FE');
    point2DTypeShiftBits = 1;

    point3DNumEvents = 0;
    point3DValid	= false(0);
    point3DTypeMask = hex2dec('FE');
    point3DTypeShiftBits = 1;
end

cellFind = @(string)(@(cellContents)(strcmp(string, cellContents)));

% Go back to the beginning of the data
fseek(fileHandle, info.beginningOfDataPointer, 'bof');

% If the file has been indexed or partially indexed, and there is a
% startPacket or startTime parameter, then jump ahead to the right place
if isfield(info, 'packetPointers')
    if startPacket > 1
        fseek(fileHandle, double(packetPointers(startPacket)), 'bof');
        packetCount = startPacket - 1;
    elseif startTime > 0
        targetPacketIndex = find(info.packetTimeStamps < startTime * 1e6, 1, 'last');
        if ~isempty(targetPacketIndex)
            fseek(fileHandle, double(packetPointers(targetPacketIndex)), 'bof');
            packetCount = targetPacketIndex - 1;
        end
    end
end

% If the file has already been indexed (PARTIAL INDEXING NOT HANDLED), and
% we are using modPacket to skip a proportion of the data, then use this
% flag to speed up the loop
modSkipping = isfield(info, 'packetPointers') && modPacket > 1;

while true % implement the exit conditions inside the loop - allows to distinguish between different types of exit
%% Headers
    % Read the header of the next packet
    packetCount = packetCount + 1;
    if modSkipping
        packetCount = ceil(packetCount / modPacket) * modPacket;
        fseek(fileHandle, double(packetPointers(packetCount)), 'bof');
    end
	header = uint8(fread(fileHandle, 28));

	if feof(fileHandle)
        packetCount = packetCount - 1;
		info.numPackets = packetCount;
		break
	end
	if length(packetTypes) < packetCount
		% Double the size of packet index arrays as necessary
		packetTypes		= [packetTypes;		ones(packetCount, 1, 'uint16') * 32768];
		packetPointers	= [packetPointers;	zeros(packetCount, 1, 'uint64')];
		packetTimeStamps	= [packetTimeStamps;	zeros(packetCount, 1, 'uint64')];
	end
	packetPointers(packetCount) = ftell(fileHandle) - 28;
	if mod(packetCount, 100) == 0
		disp(['packet: ' num2str(packetCount) '; file position: ' num2str(floor(ftell(fileHandle) / 1000000)) ' MB'])
	end
	if startPacket > packetCount || mod(packetCount, modPacket) > 0
		% Ignore this packet as its count is too low
		eventSize = typecast(header(5:8), 'int32');
		eventNumber = typecast(header(21:24), 'int32');
		fseek(fileHandle, eventNumber * eventSize, 'cof');
    elseif endPacket < packetCount
        packetCount = packetCount - 1;
		info.numPackets = packetCount;
        break
    else
		eventSize = typecast(header(5:8), 'int32');
		eventTsOffset = typecast(header(9:12), 'int32');
		eventTsOverflow = typecast(header(13:16), 'int32');
		%eventCapacity = typecast(header(17:20), 'int32');
		eventNumber = typecast(header(21:24), 'int32');
		%eventValid = typecast(header(25:28), 'int32');
		% Read the full packet
		numBytesInPacket = eventNumber * eventSize;
		packetData = uint8(fread(fileHandle, numBytesInPacket));
		% Find the first timestamp and check the timing constraints
		packetTimeStampOffset = uint64(eventTsOverflow) * uint64(2^31);
		mainTimeStamp = uint64(typecast(packetData (eventTsOffset + 1 : eventTsOffset + 4), 'int32')) + packetTimeStampOffset;
     	packetTimeStamps(packetCount) = mainTimeStamp;

        if mainTimeStamp > endTime * 1e6 && ...
                mainTimeStamp ~= hex2dec('7FFFFFFF') % This may be a timestamp reset - don't let it stop the import
            % Naively assume that the packets are all ordered correctly and finish
            packetCount = packetCount - 1;
            break
        end
        if startTime * 1e6 <= mainTimeStamp
			eventType = typecast(header(1:2), 'int16');
			packetTypes(packetCount) = eventType;

			%eventSource = typecast(data(3:4), 'int16'); % Multiple sources not handled yet

            if ~noData
    			% Handle the packet types individually:
    %% Special events
                if eventType == 0
                    if allDataTypes || any(cellfun(cellFind('special'), dataTypes))
                        % First check if the array is big enough
                        currentLength = length(specialValid);
                        if currentLength == 0
                            specialValid		= false(eventNumber, 1);
                            specialTimeStamp	= zeros(eventNumber, 1, 'uint64');
                            specialAddress		= zeros(eventNumber, 1, 'uint32');
                        else
                            while eventNumber > currentLength - specialNumEvents
                                specialValid		= [specialValid;		false(currentLength, 1)];
                                specialTimeStamp	= [specialTimeStamp;	zeros(currentLength, 1, 'uint64')];
                                specialAddress		= [specialAddress;		zeros(currentLength, 1, 'uint32')];
                                currentLength = length(specialValid);
                                %disp(['Special array resized to ' num2str(currentLength)])
                            end
                        end
                        % Iterate through the events, converting the data and
                        % populating the arrays
                        % TO DO - MATRICISE THIS COMPUTATION, FOLLOWING THE
                        % EXAMPLE IN POLARITY
                        for dataPointer = 1 : eventSize : numBytesInPacket % This points to the first byte for each event
                            specialNumEvents = specialNumEvents + 1;
                            specialValid(specialNumEvents) = mod(packetData(dataPointer), 2) == 1; %Pick off the first bit
                            specialTimeStamp(specialNumEvents) = packetTimeStampOffset + uint64(typecast(packetData(dataPointer + 4 : dataPointer + 7), 'int32'));
                            specialAddress(specialNumEvents) = uint8(bitshift(bitand(packetData(dataPointer), specialDataMask), -specialDataShiftBits));
                        end
                    end
    %% Polarity events
                elseif eventType == 1
                    if allDataTypes || any(cellfun(cellFind('polarity'), dataTypes))
                        % First check if the array is big enough
                        currentLength = length(polarityValid);
                        if currentLength == 0
                            polarityValid		= false(eventNumber, 1);
                            polarityTimeStamp	= zeros(eventNumber, 1, 'uint64');
                            polarityX			= zeros(eventNumber, 1, 'uint16');
                            polarityY			= zeros(eventNumber, 1, 'uint16');
                            polarityPolarity	= false(eventNumber, 1);
                        else
                            while eventNumber > currentLength - polarityNumEvents
                                polarityValid		= [polarityValid;		false(currentLength, 1)];
                                polarityTimeStamp	= [polarityTimeStamp;	zeros(currentLength, 1, 'uint64')];
                                polarityX			= [polarityX;			zeros(currentLength, 1, 'uint16')];
                                polarityY			= [polarityY;			zeros(currentLength, 1, 'uint16')];
                                polarityPolarity	= [polarityPolarity;	false(currentLength, 1)];
                                currentLength = length(polarityValid);
                                %disp(['Polarity array resized to ' num2str(currentLength)])
                            end
                        end
                        dataMatrix = reshape(packetData, [eventSize, eventNumber]);
                        dataTempTimeStamp = dataMatrix(5:8, :);
                        polarityTimeStamp(polarityNumEvents + (1 : eventNumber)) = packetTimeStampOffset + uint64(typecast(dataTempTimeStamp(:), 'int32'));
                        dataTempAddress = dataMatrix(1:4, :);
                        dataTempAddress = typecast(dataTempAddress(:), 'uint32');
                        polarityValid(polarityNumEvents + (1 : eventNumber)) = mod(dataTempAddress, 2) == 1; % Pick off the first bit
                        polarityPolarity(polarityNumEvents + (1 : eventNumber)) = mod(floor(dataTempAddress / 2), 2) == 1; % Pick out the second bit
                        polarityY(polarityNumEvents + (1 : eventNumber)) = uint16(bitshift(bitand(dataTempAddress, polarityYMask), -polarityYShiftBits));
                        polarityX(polarityNumEvents + (1 : eventNumber)) = uint16(bitshift(bitand(dataTempAddress, polarityXMask), -polarityXShiftBits));
                        polarityNumEvents = polarityNumEvents + eventNumber;
                    end
    %% Frames
                elseif eventType == 2
                    if allDataTypes || any(cellfun(cellFind('frame'), dataTypes))
                        % First check if the array is big enough
                        currentLength = length(frameValid);
                        if currentLength == 0
                            frameValid					= false(eventNumber, 1);
                            frameColorChannels			= zeros(eventNumber, 1, 'uint8');
                            frameColorFilter			= zeros(eventNumber, 1, 'uint8');
                            frameRoiId					= zeros(eventNumber, 1, 'uint8');
                            if simplifyFrameTimeStamps
                                frameTimeStampStart     = zeros(eventNumber, 1, 'uint64');
                                frameTimeStampEnd		= zeros(eventNumber, 1, 'uint64');
                            else
                                frameTimeStampFrameStart	= zeros(eventNumber, 1, 'uint64');
                                frameTimeStampFrameEnd		= zeros(eventNumber, 1, 'uint64');
                                frameTimeStampExposureStart = zeros(eventNumber, 1, 'uint64');
                                frameTimeStampExposureEnd	= zeros(eventNumber, 1, 'uint64');
                            end
                            frameXLength				= zeros(eventNumber, 1, 'uint16');
                            frameYLength				= zeros(eventNumber, 1, 'uint16');
                            frameXPosition				= zeros(eventNumber, 1, 'uint16');
                            frameYPosition				= zeros(eventNumber, 1, 'uint16');
                            frameSamples				= cell(eventNumber, 1);
                        else
                            while eventNumber > currentLength - frameNumEvents
                                frameValid					= [frameValid;                  false(currentLength, 1)];
                                frameColorChannels			= [frameColorChannels;			zeros(currentLength, 1, 'uint8')];
                                frameColorFilter			= [frameColorFilter;			zeros(currentLength, 1, 'uint8')];
                                frameRoiId					= [frameRoiId;					zeros(currentLength, 1, 'uint8')];
                                if simplifyFrameTimeStamps
                                    frameTimeStampStart	= [frameTimeStampStart;         zeros(currentLength, 1, 'uint64')];
                                    frameTimeStampEnd		= [frameTimeStampEnd;		zeros(currentLength, 1, 'uint64')];
                                else
                                    frameTimeStampFrameStart	= [frameTimeStampFrameStart;	zeros(currentLength, 1, 'uint64')];
                                    frameTimeStampFrameEnd		= [frameTimeStampFrameEnd;		zeros(currentLength, 1, 'uint64')];
                                    frameTimeStampExposureStart = [frameTimeStampExposureStart; zeros(currentLength, 1, 'uint64')];
                                    frameTimeStampExposureEnd	= [frameTimeStampExposureEnd;	zeros(currentLength, 1, 'uint64')];
                                end
                                frameXLength				= [frameXLength;				zeros(currentLength, 1, 'uint16')];
                                frameYLength				= [frameYLength;				zeros(currentLength, 1, 'uint16')];
                                frameXPosition				= [frameXPosition;				zeros(currentLength, 1, 'uint16')];
                                frameYPosition				= [frameYPosition;				zeros(currentLength, 1, 'uint16')];
                                frameSamples				= [frameSamples;				cell(currentLength, 1)];
                                currentLength = length(frameValid);
                                %disp(['Frame array resized to ' num2str(currentLength)])
                            end
                        end

                        % Iterate through the events, converting the data and
                        % populating the arrays
                        for dataPointer = 1 : eventSize : numBytesInPacket % This points to the first byte for each event
                            frameNumEvents = frameNumEvents + 1;
                            frameValid(frameNumEvents) = mod(packetData(dataPointer), 2) == 1; % Pick off the first bit
                            frameData = typecast(packetData(dataPointer : dataPointer + 3), 'uint32');
                            frameColorChannels(frameNumEvents) = uint16(bitshift(bitand(frameData, frameColorChannelsMask), -frameColorChannelsShiftBits));
                            frameColorFilter(frameNumEvents)	= uint16(bitshift(bitand(frameData, frameColorFilterMask),	-frameColorFilterShiftBits));
                            frameRoiId(frameNumEvents)		= uint16(bitshift(bitand(frameData, frameRoiIdMask),		-frameRoiIdShiftBits));
                            if simplifyFrameTimeStamps
                                frameTimeStampStart(frameNumEvents)		= packetTimeStampOffset + uint64(typecast(packetData(dataPointer + 12 : dataPointer + 15), 'int32'));
                                frameTimeStampEnd(frameNumEvents)		= packetTimeStampOffset + uint64(typecast(packetData(dataPointer + 16 : dataPointer + 19), 'int32'));
                            else
                                frameTimeStampFrameStart(frameNumEvents)		= packetTimeStampOffset + uint64(typecast(packetData(dataPointer + 4 : dataPointer + 7), 'int32'));
                                frameTimeStampFrameEnd(frameNumEvents)		= packetTimeStampOffset + uint64(typecast(packetData(dataPointer + 8 : dataPointer + 11), 'int32'));
                                frameTimeStampExposureStart(frameNumEvents)	= packetTimeStampOffset + uint64(typecast(packetData(dataPointer + 12 : dataPointer + 15), 'int32'));
                                frameTimeStampExposureEnd(frameNumEvents)		= packetTimeStampOffset + uint64(typecast(packetData(dataPointer + 16 : dataPointer + 19), 'int32'));
                            end
                            frameXLength(frameNumEvents)		= typecast(packetData(dataPointer + 20 : dataPointer + 21), 'uint16'); % strictly speaking these are 4-byte signed integers, but there's no way they'll be that big in practice
                            frameYLength(frameNumEvents)		= typecast(packetData(dataPointer + 24 : dataPointer + 25), 'uint16');
                            frameXPosition(frameNumEvents)	= typecast(packetData(dataPointer + 28 : dataPointer + 29), 'uint16');
                            frameYPosition(frameNumEvents)	= typecast(packetData(dataPointer + 32 : dataPointer + 33), 'uint16');
                            numSamples = int32(frameXLength(frameNumEvents)) * int32(frameYLength(frameNumEvents)) * int32(frameColorChannels(frameNumEvents)); % Conversion to int32 allows addition with 'dataPointer' below
                            % At least one recording has a file ending half way
                            % through the frame data due to a laptop dying,
                            % hence the following check
                            if length(packetData) >= dataPointer + 35 + numSamples * 2
                                sampleData = cast(typecast(packetData(dataPointer + 36 : dataPointer + 35 + numSamples * 2), 'uint16'), 'uint16');
                                frameSamples{frameNumEvents}		= reshape(sampleData, frameColorChannels(frameNumEvents), frameXLength(frameNumEvents), frameYLength(frameNumEvents));
                                if frameColorChannels(frameNumEvents) == 1
                                    frameSamples{frameNumEvents} = squeeze(frameSamples{frameNumEvents});
                                    frameSamples{frameNumEvents} = permute(frameSamples{frameNumEvents}, [2 1]);
                                else
                                    % Change the dimensions of the frame array to
                                    % the standard for matlab: column, then row,
                                    % then channel number
                                    frameSamples{frameNumEvents} = permute(frameSamples{frameNumEvents}, [3 2 1]);
                                end
                            else
                                frameSamples{frameNumEvents} = zeros(frameYLength(frameNumEvents), frameXLength(frameNumEvents), frameColorChannels(frameNumEvents));
                            end
                            % aedat3 uses left-justified 16 bit samples -
                            % for consistency with aedat2, revert to
                            % fundamental 10 bit representation
                            frameSamples{frameNumEvents} = frameSamples{frameNumEvents} / 2^6;
                        end

                    end
    %% IMU6
                elseif eventType == 3
                     if allDataTypes || any(cellfun(cellFind('imu6'), dataTypes))
                        % First check if the array is big enough
                        currentLength = length(imu6Valid);
                        if currentLength == 0
                            imu6Valid			= false(eventNumber, 1);
                            imu6TimeStamp		= zeros(eventNumber, 1, 'uint64');
                            imu6AccelX			= zeros(eventNumber, 1, 'single');
                            imu6AccelY			= zeros(eventNumber, 1, 'single');
                            imu6AccelZ			= zeros(eventNumber, 1, 'single');
                            imu6GyroX			= zeros(eventNumber, 1, 'single');
                            imu6GyroY			= zeros(eventNumber, 1, 'single');
                            imu6GyroZ			= zeros(eventNumber, 1, 'single');
                            imu6Temperature     = zeros(eventNumber, 1, 'single');
                        else
                            while eventNumber > currentLength - imu6NumEvents
                                imu6Valid			= [imu6Valid;        false(currentLength, 1)];
                                imu6TimeStamp		= [imu6TimeStamp;    zeros(currentLength, 1, 'uint64')];
                                imu6AccelX			= [imu6AccelX;       zeros(currentLength, 1, 'single')];
                                imu6AccelY			= [imu6AccelY;       zeros(currentLength, 1, 'single')];
                                imu6AccelZ			= [imu6AccelZ;       zeros(currentLength, 1, 'single')];
                                imu6GyroX			= [imu6GyroX;        zeros(currentLength, 1, 'single')];
                                imu6GyroY			= [imu6GyroY;        zeros(currentLength, 1, 'single')];
                                imu6GyroZ			= [imu6GyroZ;        zeros(currentLength, 1, 'single')];
                                imu6Temperature     = [imu6Temperature;  zeros(currentLength, 1, 'single')];
                                currentLength = length(imu6Valid);
                            end
                        end
                        % Matricise computation on a packet
                        dataMatrix = reshape(packetData, [eventSize, eventNumber]);
                        dataTempTimeStamp = dataMatrix(5:8, :);
                        imu6TimeStamp(imu6NumEvents + (1 : eventNumber)) = packetTimeStampOffset + uint64(typecast(dataTempTimeStamp(:), 'int32'));
                        % The following is overkill - only need to pick out 1
                        % byte to get to the valid flag
                        dataTempAddress = dataMatrix(1:4, :);
                        dataTempAddress = typecast(dataTempAddress(:), 'uint32');
                        imu6Valid(imu6NumEvents + (1 : eventNumber)) = mod(dataTempAddress, 2) == 1; % Pick off the first bit
                        tempDataMatrix = dataMatrix(9:12, :);
                        imu6AccelX      (imu6NumEvents + (1 : eventNumber)) = typecast(tempDataMatrix(:), 'single');
                        tempDataMatrix = dataMatrix(13:16, :);
                        imu6AccelY      (imu6NumEvents + (1 : eventNumber)) = typecast(tempDataMatrix(:), 'single');
                        tempDataMatrix = dataMatrix(17:20, :);
                        imu6AccelZ      (imu6NumEvents + (1 : eventNumber)) = typecast(tempDataMatrix(:), 'single');
                        tempDataMatrix = dataMatrix(21:24, :);
                        imu6GyroX       (imu6NumEvents + (1 : eventNumber)) = typecast(tempDataMatrix(:), 'single');
                        tempDataMatrix = dataMatrix(25:28, :);
                        imu6GyroY       (imu6NumEvents + (1 : eventNumber)) = typecast(tempDataMatrix(:), 'single');
                        tempDataMatrix = dataMatrix(29:32, :);
                        imu6GyroZ       (imu6NumEvents + (1 : eventNumber)) = typecast(tempDataMatrix(:), 'single');
                        tempDataMatrix = dataMatrix(33:36, :);
                        imu6Temperature (imu6NumEvents + (1 : eventNumber)) = typecast(tempDataMatrix(:), 'single');
                        imu6NumEvents = imu6NumEvents + eventNumber;
                    end
    %% Sample
                elseif eventType == 5
                     if allDataTypes || any(cellfun(cellFind('sample'), dataTypes))
                        %{
                        sampleValid			= bool([]);
                        sampleTimeStamp		= uint64([]);
                        sampleSampleType	= uint8([]);
                        sampleSample		= uint32([]);
                        %}
                    end
    %% Ear
                elseif eventType == 6
                     if allDataTypes || any(cellfun(cellFind('ear'), dataTypes))
                            %{
                            earValid		= bool([]);
                            earTimeStamp	= uint64([]);
                            earPosition 	= uint8([]);
                            earChannel		= uint16([]);
                            earNeuron		= uint8([]);
                            earFilter		= uint8([]);
                            %}
                    end
    %% Point1D

                elseif eventType == 8
                    if allDataTypes || any(cellfun(cellFind('point1D'), dataTypes))
                        % First check if the array is big enough
                        currentLength = length(point1DValid);
                        if currentLength == 0
                            point1DValid		= false(eventNumber, 1);
                            point1DType         = zeros(eventNumber, 1, 'uint8');
                            point1DTimeStamp	= zeros(eventNumber, 1, 'uint64');
                            point1DX            = zeros(eventNumber, 1, 'single');
                        else
                            while eventNumber > currentLength - point1DNumEvents
                                point1DValid		= [point1DValid;		false(currentLength, 1)];
                                point1DType     	= [point1DType;     	zeros(currentLength, 1, 'uint8')];
                                point1DTimeStamp	= [point1DTimeStamp;	zeros(currentLength, 1, 'uint64')];
                                point1DX            = [point1DX;            zeros(currentLength, 1, 'single')];
                                currentLength = length(point1dValid);
                            end
                        end
                        % Iterate through the events, converting the data and
                        % populating the arrays
                        for dataPointer = 1 : eventSize : numBytesInPacket % This points to the first byte for each event
                            point1DNumEvents = point1DNumEvents + 1;
                            point1DValid(point1DNumEvents) = boolean(bitand(packetData(dataPointer), 1)); %Pick off the first bit
                            point1DType(point1DNumEvents) = bitshift(bitand(packetData(dataPointer), point1DTypeMask), -point1DTypeShiftBits);
                            point1DTimeStamp(point1DNumEvents) = packetTimeStampOffset + uint64(typecast(packetData(dataPointer + 8 : dataPointer + 11), 'int32'));
                            point1DX(point1DNumEvents) = typecast(packetData(dataPointer + 4 : dataPointer + 7), 'single');
                        end
                    end
    %% Point2D
                elseif eventType == 9
                    if allDataTypes || any(cellfun(cellFind('point2D'), dataTypes))
                        % First check if the array is big enough
                        currentLength = length(point2DValid);
                        if currentLength == 0
                            point2DValid		= false(eventNumber, 1);
                            point2DType         = zeros(eventNumber, 1, 'uint8');
                            point2DTimeStamp	= zeros(eventNumber, 1, 'uint64');
                            point2DX            = zeros(eventNumber, 1, 'single');
                            point2DY            = zeros(eventNumber, 1, 'single');
                        else
                            while eventNumber > currentLength - point2DNumEvents
                                point2DValid		= [point2DValid;		false(currentLength, 1)];
                                point2DType     	= [point2DType;     	zeros(currentLength, 1, 'uint8')];
                                point2DTimeStamp	= [point2DTimeStamp;	zeros(currentLength, 1, 'uint64')];
                                point2DX            = [point2DX;            zeros(currentLength, 1, 'single')];
                                point2DY            = [point2DY;            zeros(currentLength, 1, 'single')];
                                currentLength = length(point2DValid);
                            end
                        end
                        % Iterate through the events, converting the data and
                        % populating the arrays
                        for dataPointer = 1 : eventSize : numBytesInPacket % This points to the first byte for each event
                            point2DNumEvents = point2DNumEvents + 1;
                            point2DValid(point2DNumEvents) = boolean(bitand(packetData(dataPointer), 1)); %Pick off the first bit
                            point2DType(point2DNumEvents) = bitshift(bitand(packetData(dataPointer), point2DTypeMask), -point2DTypeShiftBits);
                            point2DTimeStamp(point2DNumEvents) = packetTimeStampOffset + uint64(typecast(packetData(dataPointer + 12 : dataPointer + 15), 'int32'));
                            point2DX(point2DNumEvents) = typecast(packetData(dataPointer + 4 : dataPointer + 7), 'single');
                            point2DY(point2DNumEvents) = typecast(packetData(dataPointer + 8 : dataPointer + 11), 'single');
                        end
                    end
    %% Point3D
                elseif eventType == 10
                    if allDataTypes || any(cellfun(cellFind('point3D'), dataTypes))
                        % First check if the array is big enough
                        currentLength = length(point3DValid);
                        if currentLength == 0
                            point3DValid		= false(eventNumber, 1);
                            point3DType         = zeros(eventNumber, 1, 'uint8');
                            point3DTimeStamp	= zeros(eventNumber, 1, 'uint64');
                            point3DX            = zeros(eventNumber, 1, 'single');
                            point3DY            = zeros(eventNumber, 1, 'single');
                            point3DZ            = zeros(eventNumber, 1, 'single');
                        else
                            while eventNumber > currentLength - point3DNumEvents
                                point3DValid		= [point3DValid;		false(currentLength, 1)];
                                point3DType     	= [point3DType;     	zeros(currentLength, 1, 'uint8')];
                                point3DTimeStamp	= [point3DTimeStamp;	zeros(currentLength, 1, 'uint64')];
                                point3DX            = [point3DX;            zeros(currentLength, 1, 'single')];
                                point3DY            = [point3DY;            zeros(currentLength, 1, 'single')];
                                point3DZ            = [point3DZ;            zeros(currentLength, 1, 'single')];
                                currentLength = length(point3DValid);
                            end
                        end
                        % Iterate through the events, converting the data and
                        % populating the arrays
                        for dataPointer = 1 : eventSize : numBytesInPacket % This points to the first byte for each event
                            point3DNumEvents = point3DNumEvents + 1;
                            point3DValid(point3DNumEvents) = boolean(bitand(packetData(dataPointer), 1)); %Pick off the first bit
                            point3DType(point3DNumEvents) = bitshift(bitand(packetData(dataPointer), point3DTypeMask), -point3DTypeShiftBits);
                            point3DTimeStamp(point3DNumEvents) = packetTimeStampOffset + uint64(typecast(packetData(dataPointer + 16 : dataPointer + 19), 'int32'));
                            point3DX(point3DNumEvents) = typecast(packetData(dataPointer + 4 : dataPointer + 7), 'single');
                            point3DY(point3DNumEvents) = typecast(packetData(dataPointer + 8 : dataPointer + 11), 'single');
                            point3DZ(point3DNumEvents) = typecast(packetData(dataPointer + 12 : dataPointer + 15), 'single');
                        end
                    end
                else
                    error('Unknown event type')
                end
            end
        end
	end
	if packetCount == endPacket
		break
	end
end

%% Clip data arrays to correct size

if noData == false

    outputData = struct;

    if specialNumEvents > 0
        keepLogical = false(size(specialValid, 1), 1);
        keepLogical(1:specialNumEvents) = true;
        special.valid = specialValid(keepLogical);
        special.timeStamp = specialTimeStamp(keepLogical);
        special.address = specialAddress(keepLogical);
        outputData.special = special;
    end

    if polarityNumEvents > 0
        keepLogical = false(size(polarityValid, 1), 1);
        keepLogical(1:polarityNumEvents) = true;
        polarity.valid = polarityValid(keepLogical);
        polarity.timeStamp	= polarityTimeStamp(keepLogical);
        polarity.y			= polarityY(keepLogical);
        polarity.x			= polarityX(keepLogical);
        polarity.polarity	= polarityPolarity(keepLogical);
        outputData.polarity = polarity;
    end

    if frameNumEvents > 0
        keepLogical = false(size(frameValid, 1), 1);
        keepLogical(1:frameNumEvents) = true;
        frame.valid = frameValid(keepLogical);
        frame.roiId					= frameRoiId(keepLogical);
        frame.colorChannels			= frameColorChannels(keepLogical);
        frame.colorFilter			= frameColorFilter(keepLogical);
        if simplifyFrameTimeStamps
            frame.timeStampStart	= frameTimeStampStart(keepLogical);
            frame.timeStampEnd		= frameTimeStampEnd(keepLogical);
        else
            frame.timeStampFrameStart	= frameTimeStampFrameStart(keepLogical);
            frame.timeStampFrameEnd		= frameTimeStampFrameEnd(keepLogical);
            frame.timeStampExposureStart = frameTimeStampExposureStart(keepLogical);
            frame.timeStampExposureEnd	= frameTimeStampExposureEnd(keepLogical);
        end
        frame.samples				= frameSamples(keepLogical);
        frame.xLength				= frameXLength(keepLogical);
        frame.yLength				= frameYLength(keepLogical);
        frame.xPosition				= frameXPosition(keepLogical);
        frame.yPosition				= frameYPosition(keepLogical);
        outputData.frame = frame;
    end

    if imu6NumEvents > 0
        keepLogical = false(size(imu6Valid, 1), 1);
        keepLogical(1:imu6NumEvents) = true;
        imu6.valid = imu6Valid(keepLogical);
        imu6.timeStamp	= imu6TimeStamp(keepLogical);
        imu6.gyroX		= imu6GyroX(keepLogical);
        imu6.gyroY		= imu6GyroY(keepLogical);
        imu6.gyroZ		= imu6GyroZ(keepLogical);
        imu6.accelX		= imu6AccelX(keepLogical);
        imu6.accelY		= imu6AccelY(keepLogical);
        imu6.accelZ		= imu6AccelZ(keepLogical);
        imu6.temperature = imu6Temperature(keepLogical);
        outputData.imu6 = imu6;
    end

    if sampleNumEvents > 0
        keepLogical = false(size(sampleValid, 1), 1);
        keepLogical(1:sampleNumEvents) = true;
        sample.valid = sampleValid(keepLogical);
        sample.timeStamp	= sampleTimeStamp(keepLogical);
        sample.sampleType	= sampleSampleType(keepLogical);
        sample.sample		= sampleSample(keepLogical);
        outputData.sample = sample;
    end

    if earNumEvents > 0
        keepLogical = false(size(earValid, 1), 1);
        keepLogical(1:earNumEvents) = true;
        ear.valid = earValid(keepLogical);
        ear.timeStamp	= earTimeStamp(keepLogical);
        ear.position	= earosition(keepLogical);
        ear.channel		= earChannel(keepLogical);
        ear.neuron		= earNeuron(keepLogical);
        ear.filter		= earFilter(keepLogical);
        outputData.ear = ear;
    end

    if point1DNumEvents > 0
        keepLogical = false(size(point1DValid, 1), 1);
        keepLogical(1:point1DNumEvents) = true;
        point1D.valid       = point1DValid      (keepLogical);
        point1D.type        = point1DType       (keepLogical);
        point1D.timeStamp   = point1DTimeStamp  (keepLogical);
        point1D.x           = point1DX          (keepLogical);
        outputData.point1D = point1D;
    end

    if point2DNumEvents > 0
        keepLogical = false(size(point2DValid, 1), 1);
        keepLogical(1:point2DNumEvents) = true;
        point2D.valid       = point2DValid      (keepLogical);
        point2D.type        = point2DType       (keepLogical);
        point2D.timeStamp   = point2DTimeStamp  (keepLogical);
        point2D.x           = point2DX          (keepLogical);
        point2D.y           = point2DY          (keepLogical);
        outputData.point2D = point2D;
    end

    if point3DNumEvents > 0
        keepLogical = false(size(point3DValid, 1), 1);
        keepLogical(1:point3DNumEvents) = true;
        point3D.valid       = point3DValid      (keepLogical);
        point3D.type        = point3DType       (keepLogical);
        point3D.timeStamp   = point3DTimeStamp  (keepLogical);
        point3D.x           = point3DX          (keepLogical);
        point3D.y           = point3DY          (keepLogical);
        point3D.z           = point3DZ          (keepLogical);
        outputData.point3D = point3D;
    end
end

%% Pack packet info
info.packetTypes	= packetTypes(1 : packetCount);
info.packetPointers	= packetPointers(1 : packetCount);
info.packetTimeStamps	= packetTimeStamps(1 : packetCount);

%% Calculate data volume by type

% This calculation excludes the final packet for simplicity.
% It doesn't handle partial imports or invalid data.

packetSizes = [info.packetPointers(2 : end) - info.packetPointers(1 : end - 1) - 28; 0];
info.dataVolumeByEventType = {};
eventTypesTemp = info.packetTypes;
eventTypesTemp(eventTypesTemp == 32768) = 0;
for eventType = max(eventTypesTemp): -1 : 0 % counting down means the array is only assigned once
	info.dataVolumeByEventType(eventType + 1, 1 : 2) = [{EventTypes(eventType)} sum(packetSizes(info.packetTypes == eventType))];
end

%% Pack the data into the output structure

aedat.info = info;
if exist ('outputData', 'var')
    aedat.data = outputData;
end
% the unpacked importParams should not have been changed

%% Remove invalid events

if validOnly && noData == false
	aedat = RemoveInvalidEvents(aedat);
end

%% Add NumEvents field for each data type

if noData == false
    aedat = NumEventsByType(aedat);
end

%% Find first and last time stamps

if noData == false
    aedat = FindFirstAndLastTimeStamps(aedat);
end

disp('Import finished')

author: panpanfei


text:
function [rResult] = relm_HypothesisTest(vRatesH, vRatesN, nNumberSimulation, fMagThreshold)
%RELMTEST
%script to test two earthqauake rate hypotheses using earthquake data
%
% park1=vRatesH;
% park2=vRatesN;
% clear test null;
%     xmin(i)=park1(j,1);
%     xmax(i)=park1(j,2);
%     ymin(i)=park1(j,3);
%     ymax(i)=park1(j,4);
%     zmin(i)=park1(j,5);
%     zmax(i)=park1(j,6);
%     magmin(i)=park1(j,7);
%     magmax(i)=park1(j,8);
%     lamda1(i)=park1(j,9);
%     weight(i)=park1(j,10);
% Get the numbers of observed earthquakes per bin
vNumberQuake = vRatesH(:,11);
% Get the lower magnitude-limits per bin
vMagMin = vRatesH(:,7);
% Get the forecasted numbers of events per bin
vLambdaH = vRatesH(:,9);
vLambdaN = vRatesN(:,9);
% Get the weightings per bins
vWeightH = vRatesH(:,10);
vWeightN = vRatesN(:,10);
vWeightCombined = vWeightH .* vWeightN .* (vMagMin > fMagThreshold);

% Remove rows of matrix for which weight is zero
[nRow, nColumn] = size(vRatesH);
nNewIndex = 0;
for nCnt = 1:nRow
  if vWeightCombined(nCnt)>0
    nNewIndex = nNewIndex + 1;
    vWeight(nNewIndex) = vWeightCombined(nCnt);
    vNumberQuakeSel(nNewIndex) = vNumberQuake(nCnt);
    vLambdaHSel(nNewIndex) = vLambdaH(nCnt);
    vLambdaNSel(nNewIndex) = vLambdaN(nCnt);
    mmin(nNewIndex) = vMagMin(nCnt);
  end
end

% Get the number of events (weighted)
vNumberQuake = vWeight .* vNumberQuakeSel;
nNumberQuake = sum(vNumberQuake);

% Weight the important columns
vLambdaH = vWeight .* vLambdaHSel;
vLambdaN = vWeight .* vLambdaNSel;
% Garbage collection
clear vRatesH vRatesN vLambdaHSel vLambdaNSel vNumberQuakeSel vMagMin vWeightCombined vWeightH vWeightN;

%make a weighted magnitude-frequency plot
%
% mf=[mmin;vNumberQuake;vLambdaH;vLambdaN]';
% mfsort=sortrows(mf);
% mag=mfsort(:,1);
%
% Fobs=flip(cumsum(flip(mfsort(:,2))));
% Fth1=flip(cumsum(flip(mfsort(:,3))));
% Fth2=flip(cumsum(flip(mfsort(:,4))));
% figure%(1)
% semilogy(mag,Fobs,'r',mag,Fth1,'g',mag,Fth2,'b');
% grid;
% axis([3,8,.0001,100]);
%
%    Evaluate whether total number of quakes is consistent with H1
%
%
Nhat=sum(vLambdaH);
peq=poisspdf(nNumberQuake, Nhat); % probability of exactly Nquake
Ple=poisscdf(nNumberQuake, Nhat); % probability of less than or equal to Nquake
Pless=Ple-peq;              % probability of less than Nquake
Pmore=1-Ple;                 % probability of more than Nquake
rResult.P_H_Equal = peq;
rResult.P_H_Less = Pless;
rResult.P_H_More = Pmore;
rResult.Nhat_H = Nhat;
lamcum1=cumsum(vLambdaH)/rResult.Nhat_H;
%   Evaluate whether total number of quakes is consistent with H2
Nhat=sum(vLambdaN);
peq=poisspdf(nNumberQuake, Nhat); % probability of exactly Nquake
Ple=poisscdf(nNumberQuake, Nhat); % probability of less than or equal to Nquake
Pless=Ple-peq;              % probability of less than Nquake
Pmore=1-Ple;                 % probability of more than Nquake
rResult.P_N_Equal = peq;
rResult.P_N_Less = Pless;
rResult.P_N_More = Pmore;
rResult.Nhat_N = Nhat;
lamcum2=cumsum(vLambdaN)/rResult.Nhat_N;
%
%   simulate catalogs according to H1,
%   and evaluate likelihood scores of nsquake1 and real catalog using lamda1 and lamda2
%
try
  nsquake=simulate(nNumberQuake, vLambdaH, nNumberSimulation);
  [rResult.LLR_H, rResult.fRank11, rResult.fRank12] = Rtest(vLambdaH, vLambdaN, vNumberQuake, nsquake, vWeight);
catch
  rResult.LLR_H = nan;
  rResult.fRank11 = nan;
  rResult.fRank12 = nan;
end
%
%   simulate catalogs according to H2,
%   and evaluate likelihood scores of nsquake1 and real catalog using lamda1 and lamda2
%
try
  nsquake=simulate(nNumberQuake, vLambdaN, nNumberSimulation);
  [rResult.LLR_N, rResult.fRank21, rResult.fRank22] = Rtest(vLambdaH, vLambdaN, vNumberQuake, nsquake, vWeight);
catch
  rResult.LLR_N = nan;
  rResult.fRank21 = nan;
  rResult.fRank22 = nan;
end
%
%Plot cumulative likelihood scores for two hypotheses
%
rResult.fAlpha = sum(rResult.LLR_N > 0)/nNumberSimulation;
rResult.fBeta = sum(rResult.LLR_H < 0)/nNumberSimulation;
%index=[1:nNumberSimulation]/nNumberSimulation;
%x=[0,0];y=[0,1];
%figure_w_normalized_uicontrolunits(2);
%plot(rResult.LLR_H,index,'g',rResult.LLR_N,index,'r',x,y,'b')';
%xlabel('Likelihood ratio (Variable b/Constant b)')';
%ylabel('Fraction of cases');
%title('Green assumes variable-b hypothesis; Red assumes constant=b hypothesis');
% nNumberQuake, Nhat1,Nhat2,P1_less,P1_more,P2_less,P2_more, alpha, beta, rank11,rank12, rank21, rank22
rResult.nNumberQuake = nNumberQuake;
rResult.nNumberSimulation = nNumberSimulation;

author: CelsoReyes


text:
function candoit = canGetHessian(problem)
% Checks whether the Hessian can be computed for a problem structure.
%
% function candoit = canGetHessian(problem)
%
% Returns true if the Hessian of the cost function can be computed given
% the problem description, false otherwise.
%
% See also: canGetCost canGetDirectionalDerivative canGetGradient

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:

    candoit = isfield(problem, 'hess') || ...
              (isfield(problem, 'ehess') && canGetEuclideanGradient(problem));

    if ~candoit && ...
           (isfield(problem, 'ehess') && ~canGetEuclideanGradient(problem))
        warning('manopt:canGetHessian', ...
               ['If the Hessian is supplied as a Euclidean Hessian (ehess), ' ...
                'then the Euclidean gradient must also be supplied (egrad).']);
    end

end

author: andrewssobral

repo: lrslibrary


text:
% change in min distance from any point on the ellipse of one fly to the
% nose of another fly
function [data,units] = compute_ddell2nose(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);
  if trx(fly).nframes <= 1,
    data{i} = [];
  else
    data{i} = diff(trx(fly).dell2nose,1,2) ./ trx(fly).dt;
  end
end
units = parseunits('mm/s');

author: kristinbranson

repo: JAABA


text:
function varargout = guipas(varargin)
% GUIPAS Application M-file for guifastr.fig
%    FIG = GUIPAS launch guifastr GUI.
%    GUIPAS('callback_name', ...) invoke the named callback.

% Last Modified by GUIDE v2.5 14-Dec-2004 16:11:40

if nargin == 1  % LAUNCH GUI

	pasguifig = openfig(mfilename,'new');

	% Generate a structure of handles to pass to callbacks, and store it.
	handles = guihandles(pasguifig);
    handles.param.method='';
    handles.param.qrsevent='';
    handles.param.npc=get(handles.npc,'string');
    handles.opstatus=0;
    guidata(pasguifig, handles);

    % Intialize pop-up menu
    guipas('qrs_popup_Callback',handles.qrs_popup,[],handles,varargin{1})
    % ---------------------

    % Wait for callbacks to run and window to be dismissed:
	uiwait(pasguifig);

	if nargout > 0
		varargout{1} = pasguifig;
	end

elseif ischar(varargin{1}) % INVOKE NAMED SUBFUNCTION OR CALLBACK

	try
		if (nargout)
			[varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
		else
			feval(varargin{:}); % FEVAL switchyard
		end
	catch
		disp(lasterr);
	end

end

%| ABOUT CALLBACKS:
%| GUIDE automatically appends subfunction prototypes to this file, and
%| sets objects' callback properties to call them through the FEVAL
%| switchyard above. This comment describes that mechanism.
%|
%| Each callback subfunction declaration has the following form:
%| <SUBFUNCTION_NAME>(H, EVENTDATA, HANDLES, VARARGIN)
%|
%| The subfunction name is composed using the object's Tag and the
%| callback type separated by '_', e.g. 'slider2_Callback',
%| 'figure1_CloseRequestFcn', 'axis1_ButtondownFcn'.
%|
%| H is the callback object's handle (obtained using GCBO).
%|
%| EVENTDATA is empty, but reserved for future use.
%|
%| HANDLES is a structure containing handles of components in GUI using
%| tags as fieldnames, e.g. handles.figure1, handles.slider2. This
%| structure is created at GUI startup using GUIHANDLES and stored in
%| the figure's application data using GUIDATA. A copy of the structure
%| is passed to each callback.  You can store additional information in
%| this structure at GUI startup, and you can change the structure
%| during callbacks.  Call guidata(h, handles) after changing your
%| copy to replace the stored original so that subsequent callbacks see
%| the updates. Type "help guihandles" and "help guidata" for more
%| information.
%|
%| VARARGIN contains any extra arguments you have passed to the
%| callback. Specify the extra arguments by editing the callback
%| property in the inspector. By default, GUIDE sets the property to:
%| <MFILENAME>('<SUBFUNCTION_NAME>', gcbo, [], guidata(gcbo))
%| Add any extra arguments after the last argument, before the final
%| closing parenthesis.

% --- Executes on selection change in qrs_popup.
function varargout = qrs_popup_Callback(h, eventdata, handles, varargin)
if nargin > 3
    set(h,'string',varargin{1});
    handles.param.qrsevent=get(h,'value');
    guidata(h, handles);
else
    handles.param.qrsevent=get(h,'value');
    guidata(h, handles);
end

% --- Executes on button press in rapco_method_check.
function varargout = rapco_method_check_Callback(h, eventdata, handles, varargin)
set(handles.gmean_method_check,'value',0);
set(handles.median_method_check,'value',0);
set(handles.mean_method_check,'value',0);
set(handles.npc_text,'ForeGroundColor',[0 0 0.55]);
set(handles.npc,'Enable','on');
handles.param.method='obs';
handles.param.npc=get(handles.npc,'string');
guidata(h, handles);

% --- Executes on editing npc textbox
function varargout = npc_Callback(h, eventdata, handles, varargin)
handles.param.npc=get(handles.npc,'string');
guidata(h, handles);

% --- Executes on button press in mean_method_check.
function varargout = mean_method_check_Callback(h, eventdata, handles, varargin)
set(handles.gmean_method_check,'value',0);
set(handles.median_method_check,'value',0);
set(handles.rapco_method_check,'value',0);
set(handles.npc_text,'ForeGroundColor',[0.55 0.55 0.55]);
set(handles.npc,'Enable','off');
handles.param.method='mean';
guidata(h, handles);

% --- Executes on button press in gmean_method_check.
function varargout = gmean_method_check_Callback(h, eventdata, handles, varargin)
set(handles.mean_method_check,'value',0);
set(handles.median_method_check,'value',0);
set(handles.rapco_method_check,'value',0);
set(handles.npc_text,'ForeGroundColor',[0.55 0.55 0.55]);
set(handles.npc,'Enable','off');
handles.param.method='gmean';
guidata(h, handles);

% --- Executes on button press in median_method_check.
function varargout = median_method_check_Callback(h, eventdata, handles, varargin)
set(handles.mean_method_check,'value',0);
set(handles.gmean_method_check,'value',0);
set(handles.rapco_method_check,'value',0);
set(handles.npc_text,'ForeGroundColor',[0.55 0.55 0.55]);
set(handles.npc,'Enable','off');
handles.param.method='median';
guidata(h, handles);

% --- Executes on button press in ok_button.
function varargout = ok_button_Callback(h, eventdata, handles, varargin)
handles.param.qrsevent=get(handles.qrs_popup,'value');
if isempty(handles.param.method)
    errordlg('Please select artifact construction method.','fmrib_pas error');
    return;
end
handles.opstatus=1;
guidata(h,handles);
uiresume(handles.pasgui_fig);

% --- Executes on button press in help_button.
function varargout = help_button_Callback(h, eventdata, handles, varargin)
pophelp('pop_fmrib_pas.m');

% --- Executes on button press in cancel_button.
function varargout = cancel_button_Callback(h, eventdata, handles, varargin)
handles.opstatus=0;
guidata(h,handles);
uiresume(handles.pasgui_fig);

% --- Executes during object creation, after setting all properties.
function npc_CreateFcn(hObject, eventdata, handles)

author: PatternRecognition

repo: OpenBMI


text:
function exact = p49_exact ( )

%*****************************************************************************80
%
%% P49_EXACT returns the exact integral for problem 49.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the value of the integral.
%
  exact = 61.0 * log ( 2.0 ) + 77.0 * log ( 7.0 ) / 4.0 - 27.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = callback_active_contour(x, options)

% callback_active_contour - callback for conjugate gradient
%
%   y = callback_active_contour(x, options);
%
%   Copytight (c) 2007 Gabriel Peyre

% norm of gradient
dt = options.dt;
n = options.n;
d = options.d;
x = reshape(x,n,n);
y = divgrad( divgrad(x,options)./repmat(d,[1 1 2]),options );
if isempty(options.E)
    y = x(:) - dt * d(:) .* y(:);
else
    y = x(:) - dt * options.E(:) .* d(:) .* y(:);
end

author: alecjacobson

repo: gptoolbox


text:
classdef Compartment < handle
    properties
        count             %number of compartments
        wholeCellModelIDs %ids of compartments
        names             %names of compartments
    end

    %indices
    properties
        cytosolIndexs                   %index within compartments
        chromosomeIndexs                %index within compartments
        membraneIndexs                  %index within compartments
        terminalOrganelleCytosolIndexs  %index within compartments
        terminalOrganelleMembraneIndexs %index within compartments
        extracellularIndexs             %index within compartments
        cellularIndexs                  %index within compartments
    end

    methods
        function this = Compartment()
        end

        function initializeConstants(this, knowledgeBase, ~)
            this.count             = knowledgeBase.numCompartments;
            this.wholeCellModelIDs = {knowledgeBase.compartments.wholeCellModelID}';
            this.names             = {knowledgeBase.compartments.name}';

            this.cytosolIndexs                   = knowledgeBase.cytosolCompartmentIndexs;
            this.chromosomeIndexs                = knowledgeBase.chromosomeCompartmentIndexs;
            this.membraneIndexs                  = knowledgeBase.membraneCompartmentIndexs;
            this.terminalOrganelleCytosolIndexs  = knowledgeBase.terminalOrganelleCytosolCompartmentIndexs;
            this.terminalOrganelleMembraneIndexs = knowledgeBase.terminalOrganelleMembraneCompartmentIndexs;
            this.extracellularIndexs             = knowledgeBase.extracellularCompartmentIndexs;
            this.cellularIndexs                  = knowledgeBase.cellularCompartmentIndexs;
        end

        function index = getIndexs(this, wholeCellModelIDs)
            [~, index] = ismember(wholeCellModelIDs, this.wholeCellModelIDs);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [E,GReuss,GHill] = shearModulus(S,h,u)
% shear modulus for an compliance tensor
%
% Syntax
%
%   [GV,GR,GVRH] = shearModulus(S) % the isotropic case
%
%   E = shearModulus(S,h,u) % the anisotropic case with plane h and shear direction  u
%   E = shearModulus(S,[],u)
%   E = shearModulus(S,h,[])
%
% Input
%  C - elastic @stiffnessTensor
%  h - shear plane @vector3d
%  u - shear direction @vector3d
%
% Output
%  E - shear modulus
%  GVoigt - Voigt effective shear modulus, upper bound
%  GReuss - Reuss effective shear modulus, lower bound
%  GHill - Hill effective shear modulus
%
% Description
%
% $$E = \frac{1}{4 S_{ijkl} h_i u_j h_k u_l}$$
%
% See also
% complianceTensor/YoungsModulus complianceTensor/volumeCompressibility complianceTensor/ChristoffelTensor

if nargin == 1 % the isotropic case

  % compute stifness tensor as 6x6 matrices
  C = matrix(inv(S),'voigt');
  S = matrix(S,'voigt');

  % the Voigt upper bound
  % GV = ((C(1,1)+C(2,2)+C(3,3))-(C(1,2)+C(2,3)+C(3,1))+3*(C(4,4)+C(5,5)+C(6,6)))/15
  GVoigt = ((C(1,1,:) + C(2,2) + C(3,3,:)) ...
    - (C(1,2,:) + C(2,3,:) + C(3,1,:)) ...
    + 3 * (C(4,4,:) + C(5,5,:) + C(6,6,:))) ./ 15;

  % the Reuss lower bound
  % GR = 15/(4*(S(1,1)+S(2,2)+S(3,3))-4*(S(1,2)+S(2,3)+S(3,1))+3*(S(4,4)+S(5,5)+S(6,6)))
  GReuss = 15 ./ (4 * (S(1,1,:) + S(2,2,:) + S(3,3,:)) ...
    - 4 * (S(1,2,:) + S(2,3,:) + S(3,1,:)) ...
    + 3 * (S(4,4,:) + S(5,5,:) + S(6,6,:)));

  % Voigt Reuss Hill average
  GHill = 0.5.*(GVoigt + GReuss);
  E = GVoigt;

elseif nargin == 2 || isempty(h)

  E = S2FunHarmonicSym.quadrature(@(u) shearModulus(S,h,u),'bandwidth',4,S.CS);

elseif isempty(u)

  E = S2FunHarmonicSym.quadrature(@(u) shearModulus(S,h,u),'bandwidth',4,S.CS);

else

  % the anisotropic shear modulus
  E = 0.25./EinsteinSum(S,[-1 -2 -3 -4],h,-1,u,-2,h,-3,u,-4);

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function theTable = SQL_Add(addWhat,inputFile,forDatabase,beVocal)
% SQL_Add   Interpret a structured input file of time series, operations,
%           or master operations.
%
% By default adds the results to a linked mySQL database.
%
%---INPUTS:
% addWhat: 'mops' (for master operations), 'ops' (for operations), or 'ts'
%             (for time series)
% inputFile: the filename of the tab-delimited textfile to be read in [default
%            = INP_ts.txt or INP_ops.txt or INP_mops.txt]
%            The input file should be formatted with whitespace as a delimiter
%            between the entries to import.
% forDatabase: if true (default) write the results of interpreting the input file
%               to the default linked mySQL database for hctsa.
% beVocal: if true (default) gives user feedback on the input process.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
%% Check inputs, set defaults:
% ------------------------------------------------------------------------------

% addWhat: 'ts', 'mops', or 'ops'
if nargin < 1 || isempty(addWhat) || ~ismember(addWhat,{'ops','ts','mops'})
    error(['Error setting first input argument -- should be ''ts'' for TimeSeries ' ...
                ', ''ops'' for Operations, or ''mops'' for Master Operations']);
end

% inputFile
if nargin < 2 || isempty(inputFile)
    % Default filenames:
    switch inputFile
    case 'ts'
        inputFile = 'INP_ts.txt';
    case 'ops'
        inputFile = 'INP_ops.txt';
    case 'mops'
        inputFile = 'INP_mops.txt';
    end
end
if ~exist(inputFile,'file')
    error('Unknown file ''%s''',inputFile);
end

if nargin < 3 || isempty(forDatabase)
    forDatabase = false;
end
if nargin < 4
    beVocal = true; % Give user feedback by default:
end

% ------------------------------------------------------------------------------
% Display welcome message:
% ------------------------------------------------------------------------------
if beVocal
    fprintf(1,'Using input file: %s.\n',inputFile);
end
ticker = tic;

% ------------------------------------------------------------------------------
%% Open database connection
% ------------------------------------------------------------------------------
if forDatabase
    [dbc,databaseName] = SQL_OpenDatabase();
end

% ------------------------------------------------------------------------------
% Define strings to unify the different strands of code for time series /
% operations
% ------------------------------------------------------------------------------
switch addWhat
    case 'ts'
        % Check that the time series table exists in the database:
        theWhat = 'time series';
        theid = 'ts_id';
        thekid = 'tskw_id';
        theTable = 'TimeSeries';
        theKeywordTable = 'TimeSeriesKeywords';
        theRelTable = 'TsKeywordsRelate';
        maxL = 50000; % the longest time series length accepted in the database
    case 'ops'
        theWhat = 'operations';
        theid = 'op_id';
        thekid = 'opkw_id';
        theTable = 'Operations';
        theKeywordTable = 'OperationKeywords';
        theRelTable = 'OpKeywordsRelate';
    case 'mops'
        theWhat = 'master operations';
        theid = 'mop_id';
        theTable = 'MasterOperations';
end

% ------------------------------------------------------------------------------
% Check that the table exists in the datbase
% ------------------------------------------------------------------------------
if forDatabase
    existString = ['SHOW TABLES LIKE ''' theTable ''''];
    [output,emsg] = mysql_dbquery(dbc,existString);
    if isempty(output) % Table doesn't exist
        error(['Table %s doesn''t exist in the database %s.\n' ...
                    'Use: (1) install.m to set the database system up from scratch,\n' ...
                    '(2) SQL_CreateAllTables to create empty tables that can later be filled with custom libraries of operations,\n' ...
                    '(3) SQL_Reset to drop all tables in the database, and repopulate all tables with the default library of operations.'],...
                    theTable,databaseName);
    end
end

% ------------------------------------------------------------------------------
%% Open and read the input file
% ------------------------------------------------------------------------------

% Determine if input is a .mat file or an input text file based on filename extension:
if strcmp(inputFile(end-3:end),'.mat');
    isMatFile = true;
else
    isMatFile = false;
end

if ~isMatFile
    % Specified a input text file
    fid = fopen(inputFile);
    if (fid==-1)
        error('Could not load the specified input text file ''%s''',inputFile)
    end
    switch addWhat
    case 'ts' % Read the time series input file:
        if beVocal
            fprintf(1,['Need to format %s (Time Series input file) as: Name ' ...
                                                                'Keywords\n'],inputFile);
            fprintf(1,'Assuming no header line\n');
            fprintf(1,'Use whitespace as a delimiter and \\n for new lines...\n');
            fprintf(1,'(Be careful that no additional whitespace is in any fields...)\n');
        end
        formatSpec = '%s%s';
    case 'ops' % Read the operations input file:
        if beVocal
            fprintf(1,['Need to format %s (Operations input file) as: OperationCode ' ...
                                            'OperationName OperationKeywords\n'],inputFile);
            fprintf(1,'Assuming no header lines\n');
            fprintf(1,'Use whitespace as a delimiter and \\n for new lines...\n');
            fprintf(1,'(Be careful that no additional whitespace is in any fields...)\n');
        end
        formatSpec = '%s%s%s';
    case 'mops' % Read the master operations input file:
        if beVocal
            fprintf(1,'Need to format %s (Master Operations input file) as: MasterCode MasterLabel\n',inputFile);
            fprintf(1,'Assuming no header lines\n');
            fprintf(1,'Use whitespace as a delimiter and \\n for new lines...\n');
            fprintf(1,'(Be careful that no additional whitespace is in any fields...)\n');
        end
        formatSpec = '%s%s';
    end
    dataIn = textscan(fid,formatSpec,'CommentStyle','#','EndOfLine','\r\n','CollectOutput',true);
    fclose(fid);

    % ------------------------------------------------------------------------------
    % Show the user what's been imported:
    % ------------------------------------------------------------------------------
    dataIn = dataIn{1}; % Collect one big matrix of cells
    numItems = size(dataIn,1); % Number of items in the input file

    if numItems == 0
        error('The input file ''%s'' seems to be empty??',inputFile)
    end

    if beVocal
        fprintf(1,'Found %u %s in %s, I think. Take a look:\n\n',numItems,theWhat,inputFile);
        switch addWhat
        case 'ts'
            fprintf(1,'%s\t%s\n','-Name-','-Keywords-');
            fprint_ts = @(x)fprintf('%s\t%s\n',dataIn{x,1},dataIn{x,2});
        case 'ops'
            fprintf(1,'%s\t%s\t%s\n','-Operation Name-','-Master Label-','-Operation Keywords-');
            fprint_ops = @(x) fprintf('%s\t%s\t%s\n',dataIn{x,1},dataIn{x,2},dataIn{x,3});
        case 'mops'
            fprintf(1,'%s\t%s\n','-Master Code-','-Master Label-');
            fprint_mops = @(x) fprintf('%s\t%s\n',dataIn{x,1},dataIn{x,2});
        end

        for i = 1:min(3,numItems)
            switch addWhat
            case 'ts', fprint_ts(i);
            case 'ops', fprint_ops(i);
            case 'mops', fprint_mops(i);
            end
        end

        if numItems > 3
            fprintf(1,'..................(%u).....................\n',max(numItems-6,0));
            for i = max(numItems-2,4):numItems
                switch addWhat
                case 'ts', fprint_ts(i);
                case 'ops', fprint_ops(i);
                case 'mops', fprint_mops(i);
                end
            end
        end

        fprintf(1,'\nHow does it look? Make sure the metadata matches the headings\n');

        % Ask the question:
        if strcmp(addWhat,'ts')
            if forDatabase
                reply = input(['If we go on, we will attempt to read all time series ' ...
                            'from file and add all ' ...
                            'data to the database.\n<<<Type ''y'' to continue...>>> '],'s');
            else
                reply = input(['If we go on, we will attempt to read all time series ' ...
                            'from file and add all ' ...
                            'data to HCTSA.mat\n<<<Type ''y'' to continue...>>> '],'s');
            end
        else
            if forDatabase
                reply = input(['If we go on, we will attempt to add all ' ...
                            'data to the database.\n<<<Type ''y'' to continue...>>> '],'s');
            else
                reply = input(['If we go on, we will attempt to add all ' ...
                            'data to HCTSA.mat\n<<<Type ''y'' to continue...>>> '],'s');
            end
        end

        if ~strcmp(reply,'y')
            fprintf(1,'I didn''t think so. Come back later...\n');
            return
        end
    end

    if beVocal
        fprintf(1,'%s read.\n',inputFile);
    end

else
    % .mat file input is only allowed when importing time series
    if ~strcmp(addWhat,'ts')
        error(['.mat file input type only supported for importing time series. ' ...
                'Please specify a text file for importing master operations or operations.']);
    end

    % Load the 3 cells specifying the data:
    if ~exist(inputFile,'file')
        error('Could not load specified input .mat file: %s',inputFile);
    end
    inputData = load(inputFile,'timeSeriesData','labels','keywords');

    % Check that they're all as they should be:
    if ~isfield(inputData,'timeSeriesData') || ~isfield(inputData,'labels') || ~isfield(inputData,'keywords') ...
                || ~(iscell(inputData.timeSeriesData) || isnumeric(inputData.timeSeriesData)) ...
                || ~iscell(inputData.labels) || ~iscell(inputData.keywords)
        error(['Incorrectly formatted input file, %s\nExpecting input file to ' ...
                    'contain: ''timeSeriesData'' (cell or matrix), ' ...
                    '''labels'' (cell), and ''keywords'' (cell).'],inputFile);
    end

    % Get number of time series (as numItems):
    if iscell(inputData.timeSeriesData)
        numItems = length(inputData.timeSeriesData);
    else
        % Specified a matrix of time-series data
        numItems = size(inputData.timeSeriesData,1);
    end
    if beVocal
        fprintf(1,['We have %u time series (rows) in your data matrix, timeSeriesData ' ...
                                    '(loaded from %s).\n'],numItems,inputFile);
        if forDatabase
            fprintf(1,['Will store time-series data from matlab file in the database to an' ...
                        ' accuracy of 6 significant figures...\n']);
        end
    end

    % Check sizes match:
    if length(inputData.labels) ~= numItems || length(inputData.keywords) ~= numItems
        fprintf(1,'%u time series, %u labels, %u keywords in %s.\n',numItems,...
                    length(inputData.labels),length(inputData.keywords),inputFile);
        error('All cells in the input file %s must be the same length.',inputFile);
    end

    % Plot some to screen
    if beVocal
        plotNum = min(5,numItems); % display as many as the first 5 time series to demonstrate
        figure('color','w','WindowStyle','docked');
        for j = 1:plotNum
            subplot(plotNum,1,j);
            if iscell(inputData.timeSeriesData)
                plot(inputData.timeSeriesData{j},'k');
            else
                plot(inputData.timeSeriesData(j,:),'k');
            end
            title(sprintf('[%u/%u] %s (%s)',j,numItems,inputData.labels{j},inputData.keywords{j}),'interpreter','none')
            if j==plotNum
                xlabel('Time (samples)')
            end
            drawnow
        end

        if forDatabase
            reply = input(sprintf(['Does this look ok for the first %u time series?\nIf we continue, ' ...
                'we will attempt to add all %u time series in the input file' ...
                ' to the database.\n<<<Type ''y'' to continue...>>> '],plotNum,numItems),'s');
        else
            reply = input(sprintf(['Does this look ok for the first %u time series?\n' ...
                                    '<<<Type ''y'' to continue...>>> '],plotNum),'s');
        end
        if ~strcmp(reply,'y')
            fprintf(1,'I didn''t think so. Come back later...\n');
            return
        end
        close
    end

    % Ok, so we have inputData.timeSeriesData, inputData.labels, and inputData.keywords
end

% ------------------------------------------------------------------------------
% Construct a table of metadta for the time series / operations /
% master operations and fill a cell, toAdd, containing mySQL INSERT commands for
% each item in the input file.
% ------------------------------------------------------------------------------
if forDatabase
    % Define an inline function to add escape strings to format mySQL queries
    if forDatabase
        esc = @RA_sqlescapestring;
        % Initialize:
        toAdd = cell(numItems,1);
    end
    if beVocal
        fprintf(1,'Preparing mySQL statements to add %u %s to the database %s...',...
                    numItems,theWhat,databaseName);
    end
end

% Prepare local IDs:
ID = (1:numItems)';

switch addWhat
case 'ts' % Prepare toAdd cell for time series
    if beVocal
        figure('color','w','WindowStyle','docked');
    end
    % Record whether data was added or not (not too long and no missing values)
    wasGood = false(1,numItems);

    % Initialize:
    Name = cell(numItems,1);
    Keywords = cell(numItems,1);
    Length = zeros(numItems,1);
    Data = cell(numItems,1);
    for j = 1:numItems
        % Assign filename and keywords strings to this time series, and load it
        if isMatFile
            Name{j} = inputData.labels{j};
            Keywords{j} = inputData.keywords{j}; % Take out inverted commas from keywords lists
            if iscell(inputData.timeSeriesData)
                Data{j} = inputData.timeSeriesData{j};
                if size(Data{j},2) > size(Data{j},1)
                    fprintf(1,'Transposing time series\n');
                    Data{j} = Data{j}';
                end
                if size(Data{j},2) ~= 1
                    error('Multivariate time-series input? Each element of timeSeriesData must be univariate');
                end
            else
                Data{j} = inputData.timeSeriesData(j,:)'; % column vector
            end
        else
            Name{j} = dataIn{j,1};
            Keywords{j} = regexprep(dataIn{j,2},'\"',''); % Take out inverted commas from keywords lists
            % Read the time series from its filename:
            try
                Data{j} = dlmread(Name{j});
            catch emsg
                fprintf(1,'%s\n',emsg.message);
                error(['\nCould not read the data file for ''%s''.',...
                        ' Check that it''s in the path.'],Name{j})
            end
        end

        % Check that label is not empty:
        if isempty(Name{j})
            warning(['\n[%u/%u] This time series is assigned an empty label' ...
                        ' and will not be added...'],j,numItems)
            beep; continue
        end

        % Assign the length of the time series
        Length(j) = length(Data{j});

        % TEST 1: Is the time series longer than the maximum allowed value?
        if Length(j) > maxL
            warning(['\n[%u/%u]%s contains %u samples, this framework can efficiently ' ...
                'deal with time series up to %u samples\n',...
                '[Note that this maximum length can be modified in SQL_Add]\n',...
                'Skipping this time series...'],...
                j,numItems,Name{j},Length(j),maxL)
            beep; continue
        end
        % TEST 2: Does the time series contains any NaN of Inf values?
        if any(isnan(Data{j})) || any(~isfinite(Data{j}))
            warning(['\n[%u/%u] The time series: ''%s'' contains special values' ...
                        ' (e.g., NaN or Inf)...\nThis time series will not be added...'], ...
                        j,numItems,Name{j})
            beep; continue
        end

        % Passed both tests! Assign wasGood = true
        wasGood(j) = true;

        % If storing in a database, need to reassign the time-series data as text
        % (which will be stored as singles in the case of a .mat file data)
        % If making for a local mat file, keep the same precision as input.
        if forDatabase
            % Want time-series data as a comma-delimited string:
            if isMatFile
                xText = sprintf('%.6g,',Data{j}); % keeps 6 figures of accuracy
                xText = xText(1:end-1); % remove trailing comma
            else
                % Read in the time series from file as strings using textscan
                fid = fopen(Name{j});
                timeSeriesData_text = textscan(fid,'%s');
                fclose(fid);
                % Turn the time series into a comma-delimited string
                timeSeriesData_text = timeSeriesData_text{1}; % cell of time series values as strings
                xText = ''; % make a comma-delimited text version of the time series as xText
                for k = 1:length(timeSeriesData_text)
                    xText = [xText,',',timeSeriesData_text{k}];
                end
                xText = xText(2:end);
            end
            Data{j} = xText;
            % Prepare the data to be added to the database in an INSERT command:
            toAdd{j} = sprintf('(''%s'',''%s'',%u,''%s'')',...
                    esc(Name{j}),esc(Keywords{j}),Length{j},Data{j});
        end

        if beVocal
            % Plot the time series
            numSubplots = min(numItems,4);
            ax = subplot(numSubplots,1,mod(j-1,numSubplots)+1);
            plot(Data{j},'-k');
            ax.XLim = [1,Length(j)];
            titleText = sprintf('[%u/%u] %s (%u), keywords = %s',j,numItems,...
                            Name{j},Length(j),Keywords{j});
            title(titleText,'interpreter','none');
            fprintf(1,'\n%s --- loaded successfully.',titleText);
            pause(5e-3); % Wait 5ms to show the plotted time series!
        end
    end

    %-------------------------------------------------------------------------------
    % All data in -> TimeSeries table
    TimeSeries = table(ID,Name,Keywords,Length,Data);

    % Check for duplicates on the name field:
    if length(unique(TimeSeries.Name)) < height(TimeSeries)
        error('Input file contains duplicate names.');
    end

    % Check whether at least some passed the quality checks
    if ~any(wasGood)
        fprintf(1,'NONE of the %u time series in the input file passed quality checks.\n',...
                                length(wasGood));
        return
    end

    if beVocal
        % Tell me about it:
        if forDatabase
            textShow = ', ready to be uploaded to the mySQL database';
        else
            textShow = '';
        end
        fprintf(1,'\nAll time-series data loaded; %u/%u passed quality tests%s.\n',...
                        sum(wasGood),length(wasGood),textShow);
        if any(~wasGood)
            input(sprintf('[List %u time series that failed... (press any key)]',sum(~wasGood)),'s');
            iNoGood = find(~wasGood);
            for i = 1:length(iNoGood)
                if ~isempty(TimeSeries.Length(iNoGood(i)))
                    lengthText = sprintf(', N = %u',TimeSeries.Length(iNoGood(i)));
                else
                    lengthText = '';
                end
                fprintf(1,'*NOT UPLOADING:* [%u] %s (%s)%s\n',iNoGood(i),...
                    TimeSeries.Name{iNoGood(i)},TimeSeries.Keywords{iNoGood(i)},...
                    lengthText);
            end
            input(sprintf('[press any key to continue to add the remaining %u time series]',...
                                                        sum(wasGood)),'s');
        end
    end

case 'mops'
    % Prepare metadata for master operations
    Code = cell(numItems,1);
    Label = cell(numItems,1);
    for j = 1:numItems
        Code{j} = dataIn{j,1};
        Label{j} = dataIn{j,2};
        if forDatabase
            toAdd{j} = sprintf('(''%s'', ''%s'')',esc(Label{j}),esc(Code{j}));
        end
    end
    % Convert to table of metadata:
    MasterOperations = table(ID,Label,Code);
    if beVocal, fprintf(1,' Done.\n'); end

case 'ops'
    % Prepare metadata table for operations
    CodeString = cell(numItems,1);
    Name = cell(numItems,1);
    Keywords = cell(numItems,1);
    Label = cell(numItems,1);
    for j = 1:numItems
        CodeString{j} = dataIn{j,1};
        Name{j} = dataIn{j,2};
        Keywords{j} = dataIn{j,3};
        Label{j} = strtok(CodeString{j},'.');
        if forDatabase
            toAdd{j} = sprintf('(''%s'', ''%s'',''%s'',''%s'')',...
                    esc(Name{j}),esc(CodeString{j}),...
                    esc(Label{j}),esc(Keywords{j}));
        end
    end
    Operations = table(ID,Name,Label,Keywords,CodeString);
    if beVocal, fprintf(1,' Done.\n'); end

    % Check for duplicates in the input file:
    [uniqueOpNames,ia] = unique(Operations.Name);
    if length(uniqueOpNames) < height(Operations)
        warning(['Input file contains %u duplicate entries, which are being removed.\n' ...
                        'Inputting %u -> %u operations...'], ...
            height(Operations)-length(uniqueOpNames),...
            height(Operations),length(uniqueOpNames));
        % Check for duplicates:
        nameTable = tabulate(Operations.Name);
        isDuplicate = find([nameTable{:,2}] > 1);
        fprintf(1,'E.g., %s is a duplicate\n',nameTable{isDuplicate(1),1});
        % Only keep the unique ones:
        Operations = Operations(ia,:);
        numItems = height(Operations);
        fprintf(1,'We now have %u operations to input...\n',numItems);
        if forDatabase
            toAdd = toAdd(ia);
        end
    end
end

% ------------------------------------------------------------------------------
% If not writing to a database, then we're done
% ------------------------------------------------------------------------------
switch addWhat
case 'ts'
    theTable = TimeSeries(wasGood,:);
    % Note that IDs are according to the original input file
case 'mops'
    theTable = MasterOperations;
case 'ops'
    theTable = Operations;
end
if ~forDatabase
    fprintf(1,'Returning a table with metadata for %u %s.\n',height(theTable),theWhat);
    return
end

%-------------------------------------------------------------------------------
% From here on is for mySQL data import:
% ------------------------------------------------------------------------------
% Check for duplicates:
if beVocal
    fprintf(1,'Checking for duplicates already in the database... ');
end
switch addWhat
case 'ts'
    existingNames = mysql_dbquery(dbc,sprintf('SELECT Name FROM TimeSeries'));
    isDuplicate = ismember(theTable.Name,existingNames);
case 'ops'
    existingOperationNames = mysql_dbquery(dbc,'SELECT Name FROM Operations');
    isDuplicate = ismember(theTable.Name,existingOperationNames);
case 'mops'
    existingMasterLabels = mysql_dbquery(dbc,'SELECT MasterLabel FROM MasterOperations');
    isDuplicate = ismember(theTable.Label,existingMasterLabels);
end
if beVocal, fprintf(1,'Done.\n'); end

% ------------------------------------------------------------------------------
% Tell the user about duplicates
% ------------------------------------------------------------------------------
if all(isDuplicate)
    fprintf(1,'All %u %s from %s already exist in %s---no new %s to add!\n',...
                    numItems,theWhat,inputFile,databaseName,theWhat);
    return
elseif sum(isDuplicate) > 0
    if beVocal
        fprintf(1,'I found %u duplicate %s already in the database %s!\n',...
                        sum(isDuplicate),theWhat,databaseName);
        fprintf(1,'There are %u new %s to add to %s.\n',...
                        sum(~isDuplicate),theWhat,databaseName);
    end
end

% ------------------------------------------------------------------------------
% Incorporate isDuplicate with wasGood
% ------------------------------------------------------------------------------
if strcmp(addWhat,'ts')
    isBad = (isDuplicate | ~wasGood);
    if all(isBad)
        fprintf(1,['All input time series either already exist in the database, ' ...
                    'or did not pass the required quality checks. Exiting.']);
        return
    end
else
    isBad = isDuplicate;
end

% ------------------------------------------------------------------------------
%% Select the maximum id already in the table
% ------------------------------------------------------------------------------
maxId = mysql_dbquery(dbc,sprintf('SELECT MAX(%s) FROM %s',theid,theTable));
if isempty(maxId) || isempty(maxId{1}) || isnan(maxId{1})
    % No time series exist in the database yet
    maxId = 0;
else
    maxId = maxId{1}; % the maximum id -- the new items will have ids greater than this
end

% ------------------------------------------------------------------------------
%% Assemble and execute the INSERT queries
% ------------------------------------------------------------------------------
fprintf('Adding %u new %s to the %s table in %s...',...
                    sum(~isBad),theWhat,theTable,databaseName)
switch addWhat
case 'ts' % Add time series to the TimeSeries table, X at a time
          % (so as not to exceed the max_allowed_packet when transmitting the
          % time-series data).
          % An appropriate chunk size will depend on the length of time series
          % in general...
          % A kind of big problem if you add some here and then fail, because
          % then there will be inconsistencies in the Results table.
          % Safer would be to initialize results as you go in case you get a
          % fail at some chunk number... :/
    maxLength = max(TimeSeries.Length);
    % 4 at a time works for 20,000 sample time series using default mySQL settings.
    % Calibrate to this:
    numPerChunk = max([floor(20000/maxLength*4),1]); % no fewer than 1 per chunk
    numPerChunk = min([numPerChunk,500]); % no more than 500 per chunk
    SQL_AddChunked(dbc,['INSERT INTO TimeSeries (Name, Keywords, Length, ' ...
                                    'Data) VALUES'],toAdd(~isBad),numPerChunk);
case 'ops' % Add operations to the Operations table 500 at a time
    SQL_AddChunked(dbc,['INSERT INTO Operations (Name, Code, MasterLabel, ' ...
                                    'Keywords) VALUES'],toAdd(~isBad),500);
case 'mops' % Add master operations to the MasterOperations table 500 at a time
    SQL_AddChunked(dbc,['INSERT INTO MasterOperations (MasterLabel, ' ...
                                'MasterCode) VALUES'],toAdd(~isBad),500);
end
fprintf(1,' Done.\n');

% ------------------------------------------------------------------------------
% Add new entries to the Results table
% ------------------------------------------------------------------------------
if ~strcmp(addWhat,'mops')
    resultsTic = tic;
    fprintf(1,'Updating the Results table in %s...(This could take a while)...',databaseName);
    if strcmp(addWhat,'ts') && height(TimeSeries) > 1000
        fprintf(1,'\n(e.g., ~10 hours for 20,000 time series with 8,000 operations---please be patient!)...')
    end
    % The following turns out to be quite slow when adding to an existing database.
    % Maybe it does the cross join before filtering on ID. So if you have a large
    % database already, it takes ages to do the cross join, only to filter out most
    % of it anyway...?
    switch addWhat
    case 'ts'
        [~,emsg] = mysql_dbexecute(dbc,sprintf(['INSERT INTO Results (ts_id,op_id) SELECT t.ts_id,o.op_id ' ...
                    'FROM TimeSeries t CROSS JOIN Operations o ON t.ts_id > %u ORDER BY t.ts_id, o.op_id'],maxId));
    case 'ops'
        [~,emsg] = mysql_dbexecute(dbc,sprintf(['INSERT INTO Results (ts_id,op_id) SELECT t.ts_id,o.op_id ' ...
                    'FROM TimeSeries t CROSS JOIN Operations o ON o.op_id > %u ORDER BY t.ts_id, o.op_id'],maxId));
    end
    if ~isempty(emsg)
        fprintf(1,' Error. This is really really not good.\n%s',emsg);
        keyboard
    else
        fprintf(1,' initialized in %s!\n',BF_TheTime(toc(resultsTic)));
    end
end

if ~strcmp(addWhat,'mops')
    % ------------------------------------------------------------------------------
    % Update the keywords table
    % ------------------------------------------------------------------------------
    fprintf(1,'Updating the %s table in %s...',theKeywordTable,databaseName);

    % First find unique keywords from new time series by splitting against commas
    switch addWhat
    case 'ts'
        kws = TimeSeries.Keywords(~isBad);
    case 'ops'
        kws = Operations.Keywords(~isBad);
    end

    % Split into each individual keyword:
    kwSplit = regexp(kws,',','split','ignorecase');
    ukws = unique(horzcat(kwSplit{:}));
    nkw = length(ukws); % The number of unique keywords in the new set of time series
    if beVocal
        fprintf(1,'\nI found %u unique keywords in the %u new %s in %s...',...
                        nkw,sum(~isBad),theWhat,inputFile);
    end

    % How many overlap with existing keywords??:
    allkws = mysql_dbquery(dbc,sprintf('SELECT Keyword FROM %s',theKeywordTable));
    if ~isempty(allkws) % the table may be empty, in which case all keywords will be new
        isNew = ~ismember(ukws,allkws);
        % cellfun(@(x)~isempty(x),regexp(ukws,allkws,'ignorecase')); % ignore case for keywords
    else
        isNew = ones(nkw,1); % All are new
    end

    if sum(isNew) > 0
        if beVocal
            fprintf(1,['\nIt turns out that %u keywords are completely new and will be added ' ...
                        'to the %s table in %s...'],sum(isNew),theKeywordTable,databaseName);
        end
        % Add the new keywords to the Keywords table
        insertString = sprintf('INSERT INTO %s (Keyword,NumOccur) VALUES',theKeywordTable);
        fisNew = find(isNew); % Indicies of new keywords
        toAdd = cell(sum(isNew),1);
        for k = 1:sum(isNew);
            toAdd{k} = sprintf('(''%s'',0)',ukws{fisNew(k)});
        end
        SQL_AddChunked(dbc,insertString,toAdd,100);
        fprintf(1,' Added %u new keywords!\n',sum(isNew));
    else
        if beVocal
            fprintf(1,['\nIt turns out that all new keywords already exist in ' ...
                    'the %s table in %s -- there are no new keywords to add\n'],...
                    theKeywordTable,databaseName);
        end
    end

    % ------------------------------------------------------------------------------
    %% Fill new keyword relationships
    % ------------------------------------------------------------------------------
    fprintf(1,'Writing new keyword relationships to the %s table in %s...', ...
                                            theRelTable,databaseName);

    % Try doing it from scratch...:
    switch addWhat
    case 'ts'
        allNames = BF_cat(TimeSeries.Name(~isBad),',','''');
    case 'ops'
        allNames = BF_cat(Operations.Name(~isBad),',','''');
    end

    % ourids: ids matching time series or operations
    ourids = mysql_dbquery(dbc,sprintf('SELECT %s FROM %s WHERE Name IN (%s)',theid, ...
                                                        theTable,allNames));
    ourids = vertcat(ourids{:});

    % ourkids: kw_ids matching keywords
    ourkids = mysql_dbquery(dbc,sprintf('SELECT %s FROM %s WHERE Keyword IN (%s)', ...
                                            thekid,theKeywordTable,BF_cat(ukws,',','''')));
    ourkids = vertcat(ourkids{:});
    addCell = {};
    for i = 1:length(kwSplit)
        for j = 1:length(kwSplit{i})
            try
                addCell{end+1} = sprintf('(%u,%u)',ourids(i),ourkids(strcmp(kwSplit{i}{j},ukws)));
            catch
                keyboard
            end
        end
    end
    % Add them all in 1000-sized chunks
    SQL_AddChunked(dbc,sprintf('INSERT INTO %s (%s,%s) VALUES',theRelTable,theid,thekid),addCell,1000);

    % Update Nmatches in the keywords table
    fprintf(1,' Done.\nNow computing all match counts for all keywords...\n');
    SQL_FlushKeywords(addWhat);
end

% ------------------------------------------------------------------------------
%% Update links between operations and master operations
% ------------------------------------------------------------------------------
if ismember(addWhat,{'mops','ops'}) % there may be new links
    % Add mop_ids to Operations table
    fprintf(1,'Evaluating links between operations and master operations...'); tic
    updateString = ['UPDATE Operations AS o SET o.mop_id = (SELECT mop_id FROM MasterOperations AS m ' ...
                        'WHERE m.MasterLabel = o.MasterLabel) WHERE mop_id IS NULL'];
    [~,emsg] = mysql_dbexecute(dbc,updateString);

    if isempty(emsg)
        fprintf(' Done.\n');
    else
        error('\nOops! Error finding links between Operations and MasterOperations:\n%s\n',emsg);
    end

    updateString = sprintf(['UPDATE MasterOperations AS m SET NPointTo = ' ...
                    '(SELECT COUNT(o.mop_id) FROM Operations AS o WHERE m.mop_id = o.mop_id)']);
    [~,emsg] = mysql_dbexecute(dbc,updateString);
    if ~isempty(emsg)
        error('Error counting NPointTo operations for mop_id = %u\n%s\n',M_ids(k),emsg);
    end

end

% ------------------------------------------------------------------------------
%% Close database
% ------------------------------------------------------------------------------
SQL_CloseDatabase(dbc)

% ------------------------------------------------------------------------------
%% Tell the user all about it
% ------------------------------------------------------------------------------
fprintf('All tasks completed in %s.\nRead %s then added %u %s into %s.\n', ...
            BF_TheTime(toc(ticker)),inputFile,sum(~isBad),theWhat,databaseName);

if strcmp(addWhat,'ts') && ~isMatFile
    fprintf(1,['**The imported data are now in %s and no longer ' ...
                    'need to be in the Matlab path**\n'],databaseName);
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function STEPWISE = stepwise_tor(dat,y,varargin)
% :Usage:
% ::
%
%     STEPWISE = stepwise_tor(dat, y, [pred. names], [alpha])
%
% Stepwise regression using Matlab with a couple of extras:
%
% Print omnibus F-values for stepwise regression
%
% get adjusted R-squared
%
% save output structure
%
% ..
%    tor wager
% ..

alph = .05;
if length(varargin) > 1
    alph = varargin{2};
end

if length(varargin) > 0
    nms = varargin{1};
else
    for i = 1:size(dat,2), nms{i} = ['V ' num2str(i)]; end
end

    % Add logistic regression for categorical 1/0 DVs!

    [STEPWISE.b,STEPWISE.se,STEPWISE.pval, ...
    STEPWISE.inmodel,stats] = ...
    stepwisefit(dat,y,'penter',alph,'display','off');

    t = stats.TSTAT;
    STEPWISE.t = t;

    S = STEPWISE;
    tabled = [S.b S.se t S.inmodel' S.pval]; tabnames = {'Beta' 'Std. Err.' 't-value' 'In Model' 'p-value'};
    print_matrix(tabled, tabnames,nms);

    N = stats.dfe+stats.df0;
    r2 = (stats.SStotal-stats.SSresid) ./ stats.SStotal;
    adjr2 = 1 - (1-r2)*((N - 1) ./ (N - stats.df0 - 1));

    fprintf(1,'\nOmnibus F(%3.0f,%3.0f) = %3.2f, RMSE = %3.2f, p = %3.6f, Adj R^2 = %3.2f\n', ...
    stats.df0,stats.dfe,stats.fstat, ...
    stats.rmse,stats.pval,adjr2);
    stats.adjr2 = adjr2;
    STEPWISE.stats = stats;

 return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore