text:
function [sFile, ChannelMat] = in_fopen_tdt(DataFile)
% IN_FOPEN_TDT Open recordings saved in the Tucker Davis Technologies format.
%
% The importer needs the folder that the files are in. K Nasiotis selected one type
% of files to work as the "raw file" - (.Tbk)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Author: Konstantinos Nasiotis 2019-2020

% Not available in the compiled version
if bst_iscompiled()
    error('Reading TDT files is not available in the compiled version of Brainstorm.');
end
% Install/load TDT-SDK library
[isInstalled, errMsg] = bst_plugin('Install', 'tdt-sdk');
if ~isInstalled
    error(errMsg);
end

 %% ===== GET FILES =====
% Get base dataset folder
[DataFolder, FileName] = bst_fileparts(DataFile);

hdr.BaseFolder = DataFolder;

 %% ===== FILE COMMENT =====
% Comment: BaseFolder
Comment = DataFolder;

%% ===== READ DATA HEADERS =====
bst_progress('text', 'TDT: Reading headers...');

% Load one second segment to see what type of signals exist in this dataset
% Use as general sampling rate the rate of the HIGHEST sampled signal
% The signals that have a lower sampling rate will be interpolated to match
% the general sampling rate

headers = TDTbin2mat(DataFolder, 'HEADERS', 1);

data = TDTbin2mat(DataFolder, 'T1', 0, 'T2', 1); % 1 second segment
all_streams = fieldnames(data.streams);

%% Get rid of streams that are pNe and SynC

all_streams = all_streams(~ismember(all_streams,{'pNe1','pNe2','pNe3','pNe4','SynC'}));

% The sampling rates present are the weirdest numbers I have ever seen:
% e.g. Fs = 3051.7578125 Hz !!!
% Those numbers create problems when loading segments of data.
% The segment loading is in TimeBounds, not SampleBounds that makes it even
% worse with those sampling rates
stream_info = struct;

LFP_label_exists = 0;

ii = 1;
for iStream = 1:length(all_streams)
    stream_info(iStream).label          = all_streams{iStream};
    stream_info(iStream).fs             = data.streams.(all_streams{iStream}).fs;
    stream_info(iStream).total_channels = size(data.streams.(all_streams{iStream}).data,1);
    stream_info(iStream).channelIndices = ii:ii+size(data.streams.(all_streams{iStream}).data,1)-1;

    ii = ii + size(data.streams.(all_streams{iStream}).data,1);

    % Brainstorm needs a single sampling rate. I assign the one on the LFPs

    theLFPlabel = 'LFP';
    if strfind(all_streams{iStream},theLFPlabel)
        data_new = TDTbin2mat(DataFolder, 'STORE', all_streams{iStream},'T1', 0, 'T2', 1); % 1 second segment

        general_sampling_rate = data_new.streams.(all_streams{iStream}).fs;
        LFP_label_exists = 1;
    end
    popup_labels{iStream} = [all_streams{iStream} ' - ' num2str(stream_info(iStream).fs) ' Hz'];
end

if ~LFP_label_exists
    [indx,tf] = listdlg('PromptString',{'Select the label of the electrophysiological signal that is present in this dataset',...
    'If more streams are present, they will be resampled to match the Fs of this stream.',''},...
    'SelectionMode','single','ListString', popup_labels);

    data_new = TDTbin2mat(DataFolder, 'STORE', all_streams{indx},'T1', 0, 'T2', 1); % 1 second segment

    general_sampling_rate = data_new.streams.(all_streams{indx}).fs;
    LFP_label_exists = 1;

    if isempty(data_new.streams)
        error('The selected stream is empty');
    end

    if isempty(indx)
        error('No stream was selected');
    end
end

nChannels = sum([stream_info.total_channels]);

 %% ===== CREATE BRAINSTORM SFILE STRUCTURE =====
% Initialize returned file structure
sFile = db_template('sfile');

% Add information read from header
sFile.prop.sfreq   =  general_sampling_rate;
sFile.byteorder    = 'l';
sFile.filename     = DataFolder;
sFile.format       = 'EEG-TDT';
sFile.device       = 'Tucker Davis Technologies';
sFile.header.tdt   = headers;
sFile.comment      = FileName;
sFile.prop.times   = [0, headers.stopTime-headers.startTime];
sFile.prop.nAvg    = 1;
% No info on bad channels
sFile.channelflag  = ones(nChannels, 1);

sFile.header.stream_info = stream_info;

%% ===== CREATE EMPTY CHANNEL FILE =====
ChannelMat = db_template('channelmat');
ChannelMat.Comment = 'TDT channels';
ChannelMat.Channel = repmat(db_template('channeldesc'), [1, nChannels]);

ii = 0;
for iStream = 1:length(all_streams)
     for iChannel = 1:stream_info(iStream).total_channels
         ii = ii+1;
         if ~(stream_info(iStream).total_channels==1)
            ChannelMat.Channel(ii).Name = [all_streams{iStream} '_' num2str(iChannel)];
         else
             ChannelMat.Channel(ii).Name= [all_streams{iStream}];
         end
         ChannelMat.Channel(ii).Loc     = [0;0;0];

         ChannelMat.Channel(ii).Group   = all_streams{iStream};

         if strfind(stream_info(iStream).label, theLFPlabel)
            ChannelMat.Channel(ii).Type = 'EEG'; % Not all are EEGs - NOT SURE WHAT TO PUT HERE - AS A STARTING POINT, I PUT WHATEVER IS ONLY ONE CHANNEL SET IT AS Misc
         else
            ChannelMat.Channel(ii).Type = all_streams{iStream};
         end
         ChannelMat.Channel(ii).Orient  = [];
         ChannelMat.Channel(ii).Weight  = 1;
         ChannelMat.Channel(ii).Comment = [];
     end
end

%% Check for acquisition events

bst_progress('text', 'TDT: Collecting acquisition events...');

disp('Getting Acquisition System events')
NO_data = TDTbin2mat(DataFolder, 'TYPE', 2); % Just load epocs / events

are_there_events = ~isempty(NO_data.epocs);

if are_there_events

    all_event_Labels = fieldnames(NO_data.epocs);

    iindex = 0;

    for iEvent = 1:length(all_event_Labels)

        if sum(ismember({'Tick','Swep','Swe+'},NO_data.epocs.(all_event_Labels{iEvent}).name))~=0
            iindex = iindex + 1;

            events(iindex).label      = NO_data.epocs.(all_event_Labels{iEvent}).name;
            events(iindex).color      = rand(1,3);
            events(iindex).epochs     = ones(1,length(NO_data.epocs.(all_event_Labels{iEvent}).onset))  ;
            events(iindex).times      = NO_data.epocs.(all_event_Labels{iEvent}).onset';
            events(iindex).reactTimes = [];
            events(iindex).select     = 1;
            events(iindex).channels   = [];
            events(iindex).notes      = [];
        else
            conditions_in_event = unique(NO_data.epocs.(all_event_Labels{iEvent}).data);

            for iCondition = 1:length(conditions_in_event)

                selected_Events_for_condition = find(NO_data.epocs.(all_event_Labels{iEvent}).data == conditions_in_event(iCondition));

                iindex = iindex+1;

                if strcmp(all_event_Labels{iEvent}, 'Note')
                    events(iindex).label  = NO_data.epocs.Note.notes{iCondition};
                else
                    events(iindex).label  = [NO_data.epocs.(all_event_Labels{iEvent}).name num2str(conditions_in_event(iCondition))];
                end
                events(iindex).color      = rand(1,3);
                events(iindex).epochs     = ones(1,length(selected_Events_for_condition))  ;
                events(iindex).times      = NO_data.epocs.(all_event_Labels{iEvent}).onset(selected_Events_for_condition)';
                events(iindex).reactTimes = [];
                events(iindex).select     = 1;
                events(iindex).channels   = [];
                events(iindex).notes      = [];
            end
        end
    end
end

%% Check for spike events

check_for_spikes = 1;

if check_for_spikes
    bst_progress('text', 'TDT: Collecting spiking events...');
    NO_data = TDTbin2mat(DataFolder, 'TYPE', 3); % Just load spikes
    are_there_spikes = ~isempty(NO_data.snips);
else
    are_there_spikes = 0;
end

%%%%%%%%
disp('***************************************************')
disp('CHECK THE SPIKES. THEY ARE ONLY ASSIGNED ON RIG TWO')
disp('***************************************************')
%%%%%%%%

if are_there_spikes

    if  ~exist ('events','var')
        events = struct;
        last_event_index = 0;
    else
        last_event_index = length(events);
    end

    all_spike_event_Labels = fieldnames(NO_data.snips);

    for iSpikeDetectedField = 1:length(all_spike_event_Labels)
%         channels_are_EEG = find(strcmp({ChannelMat.Channel.Type}, 'EEG'));

        try
            channels_are_EEG_on_selected_RIG = find(strcmp({ChannelMat.Channel.Type}, 'EEG') & strcmp({ChannelMat.Channel.Group}, ['LFP' num2str(iSpikeDetectedField)]));
        catch
            error('There is an assumption here that LFP1 corresponds to eNe1, LFP2 to eNe2 and so on. If that''s not the case on your system please contact the forum')
        end

        for iChannel = 1:length(channels_are_EEG_on_selected_RIG)

            NeuronIDs = unique(NO_data.snips.(all_spike_event_Labels{iSpikeDetectedField}).sortcode(find(NO_data.snips.(all_spike_event_Labels{iSpikeDetectedField}).chan == iChannel)));

            if length(NeuronIDs)>1 && sum(ismember(NeuronIDs,0))~=0
                warning('There are Sorted AND Unsorted Spikes in this Dataset - Probably the selection for Online sorting was made mid-Recording - THE UNSORTED SPIKES WILL BE IGNORED')
                NeuronIDs = NeuronIDs(~ismember(NeuronIDs,[0,31]));
            else
                NeuronIDs = NeuronIDs(~ismember(NeuronIDs, 31));
            end

            for iNeuron = 1:length(NeuronIDs)

                if NeuronIDs(iNeuron) ~= 31 && (length(NeuronIDs) > 1 && NeuronIDs(iNeuron) ~= 0) % MARYSE MENTIONED THAT .DATA = 31 INDICATES ARTIFACTS - NOISE

                    last_event_index = last_event_index + 1;

                    if length(NeuronIDs) == 1
                        SpikesOfThatNeuronOnChannel_Indices = find(NO_data.snips.(all_spike_event_Labels{iSpikeDetectedField}).chan == iChannel & NO_data.snips.(all_spike_event_Labels{iSpikeDetectedField}).sortcode == 0); % Unsorted
                        events(last_event_index).label = ['Spikes Channel ' ChannelMat.Channel(channels_are_EEG_on_selected_RIG(iChannel)).Name];
                    else
                        SpikesOfThatNeuronOnChannel_Indices = find(NO_data.snips.(all_spike_event_Labels{iSpikeDetectedField}).chan == iChannel & NO_data.snips.(all_spike_event_Labels{iSpikeDetectedField}).sortcode == NeuronIDs(iNeuron)); % Sorted
                        events(last_event_index).label = ['Spikes Channel ' ChannelMat.Channel(channels_are_EEG_on_selected_RIG(iChannel)).Name ' |' num2str(NeuronIDs(iNeuron)) '|'];
                    end

                    events(last_event_index).color      = rand(1,3);
                    events(last_event_index).epochs     = ones(1,length(SpikesOfThatNeuronOnChannel_Indices));
                    events(last_event_index).times      = NO_data.snips.(all_spike_event_Labels{iSpikeDetectedField}).ts(SpikesOfThatNeuronOnChannel_Indices)';
                    events(last_event_index).reactTimes = [];
                    events(last_event_index).select     = 1;
                    events(last_event_index).channels   = repmat({{ChannelMat.Channel(channels_are_EEG_on_selected_RIG(iChannel)).Name}}, 1, size(events(last_event_index).times, 2));
                    events(last_event_index).notes      = [];
                end
            end
        end
    end
end

% Import this list
sFile = import_events(sFile, [], events);

end

author: brainstorm-tools


text:
function [headmodel] = ft_datatype_headmodel(headmodel, varargin)

% FT_DATATYPE_HEADMODEL describes the FieldTrip MATLAB structure for a volume
% conduction model of the head that can be used for forward computations of
% the EEG potentials or the MEG fields. The volume conduction model represents
% the geometrical and the conductive properties of the head. These determine
% how the secondary (or impressed) currents flow and how these contribute to
% the model potential or field.
%
% A large number of forward solutions for the EEG and MEG are supported
% in FieldTrip, each with its own specification of the MATLAB structure that
% describes the volume conduction model of th ehead. It would be difficult to
% list all the possibilities here. One common feature is that the volume
% conduction model should specify its type, and that preferably it should
% specify the geometrical units in which it is expressed (e.g. mm, cm or m).
%
% An example of an EEG volume conduction model with 4 concentric spheres is:
%
% headmodel =
%        r: [86 88 94 100]
%        c: [0.33 1.00 0.042 0.33]
%        o: [0 0 0]
%     type: 'concentricspheres'
%     unit: 'mm'
%
% An example of an MEG volume conduction model with a single sphere fitted to
% the scalp with its center 4 cm above the line connecting the ears is:
%
% headmodel =
%        r: [12]
%        o: [0 0 4]
%     type: 'singlesphere'
%     unit: 'cm'
%
% For each of the methods XXX for the volume conduction model, a corresponding
% function FT_HEADMODEL_XXX exists that contains all specific details and
% references to literature that describes the implementation.
%
% Required fields:
%   - type
%
% Optional fields:
%   - unit
%
% Deprecated fields:
%   - inner_skull_surface, source_surface, skin_surface, source, skin
%
% Obsoleted fields:
%   - <none specified>
%
% Revision history:
%
% (2015/latest) Use the field name "pos" instead of "pnt" for vertex positions.
%
% (2014) All numeric values are represented in double precision.
%
% (2013) Always use the field "cond" for conductivity.
%
% (2012) Use consistent names for the volume conductor type in the structure, the
% documentation and for the actual implementation, e.g. bem_openmeeg -> openmeeg,
% fem_simbio -> simbio, concentric -> concentricspheres. Deprecated the fields
% that indicate the index of the innermost and outermost surfaces.
%
% See also FT_PREPARE_HEADMODEL, FT_DATATYPE, FT_DATATYPE_COMP, FT_DATATYPE_DIP,
% FT_DATATYPE_FREQ, FT_DATATYPE_MVAR, FT_DATATYPE_RAW, FT_DATATYPE_SOURCE,
% FT_DATATYPE_SPIKE, FT_DATATYPE_TIMELOCK, FT_DATATYPE_VOLUME

% Copyright (C) 2011-2012, Cristiano Micheli, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% get the optional input arguments, which should be specified as key-value pairs
version = ft_getopt(varargin, 'version', 'latest');

if strcmp(version, 'latest')
  version = '2015';
end

if isempty(headmodel)
  return;
end

if iscell(headmodel)
  % this might represent combined EEG, ECoG and/or MEG
  for i=1:numel(headmodel)
    % call recursively
    headmodel{i} = ft_datatype_headmodel(headmodel{i}, varargin{:});
  end
  return
end

switch version

  case '2015'
    % first make it consistent with the 2014 version
    headmodel = ft_datatype_headmodel(headmodel, 'version', '2013');

    % rename pnt into pos
    headmodel = fixpos(headmodel);

  case '2014'
    % first make it consistent with the 2013 version
    headmodel = ft_datatype_headmodel(headmodel, 'version', '2013');

    % ensure that all numbers are represented in double precision
    headmodel = ft_struct2double(headmodel);

  case '2013'
    % first make it consistent with the 2012 version
    headmodel = ft_datatype_headmodel(headmodel, 'version', '2012');

    % then rename (if necessary the c into cond
    if isfield(headmodel, 'c') && ~isfield(headmodel, 'cond')
      headmodel.cond = headmodel.c;
      headmodel = rmfield(headmodel, 'c');
    elseif isfield(headmodel, 'cond') && isfield(headmodel, 'c') && isequal(headmodel.cond, headmodel.c)
      headmodel = rmfield(headmodel, 'c');
    elseif isfield(headmodel, 'cond') && isfield(headmodel, 'c') && ~isequal(headmodel.cond, headmodel.c)
      ft_error('inconsistent specification of conductive properties for %s model', headmodel.type);
    end

  case '2012'
    % the following will be determined on the fly in ft_prepare_vol_sens
    if isfield(headmodel, 'skin_surface'),        headmodel = rmfield(headmodel, 'skin_surface');        end
    if isfield(headmodel, 'source_surface'),      headmodel = rmfield(headmodel, 'source_surface');      end
    if isfield(headmodel, 'inner_skull_surface'), headmodel = rmfield(headmodel, 'inner_skull_surface'); end
    if isfield(headmodel, 'skin'),                headmodel = rmfield(headmodel, 'skin');                end
    if isfield(headmodel, 'source'),              headmodel = rmfield(headmodel, 'source');              end

    % ensure a consistent naming of the volume conduction model types
    % these should match with the FT_HEADMODEL_XXX functions
    if isfield(headmodel, 'type')
      if strcmp(headmodel.type, 'concentric')
        headmodel.type = 'concentricspheres';
      elseif strcmp(headmodel.type, 'nolte')
        headmodel.type = 'singleshell';
      elseif strcmp(headmodel.type, 'multisphere')
        headmodel.type = 'localspheres';
      elseif strcmp(headmodel.type, 'bem_cp')
        headmodel.type = 'bemcp';
      elseif strcmp(headmodel.type, 'bem_dipoli')
        headmodel.type = 'dipoli';
      elseif strcmp(headmodel.type, 'bem_asa')
        headmodel.type = 'asa';
      elseif strcmp(headmodel.type, 'bem_openmeeg')
        headmodel.type = 'openmeeg';
      elseif strcmp(headmodel.type, 'fem_simbio')
        headmodel.type = 'simbio';
      elseif strcmp(headmodel.type, 'fdm_fns')
        headmodel.type = 'fns';
      elseif strcmp(headmodel.type, 'bem')
        ft_error('not able to convert the original ''bem'' volume type, try using headmodel.type=''dipoli''');
      elseif strcmp(headmodel.type, 'avo')
        ft_error('this format is not supported anymore');
      end
    end

    if isfield(headmodel, 'sens')
      % this applies to type=interpolate, ensure that the sensor description is up to date
      headmodel.sens = ft_datatype_sens(headmodel.sens);
    end

    if isfield(headmodel, 'type') && any(strcmp(headmodel.type, {'concentricspheres', 'singlesphere'}))
      if isfield(headmodel, 'cond') && ~isfield(headmodel, 'c')
        headmodel.c = headmodel.cond;
        headmodel = rmfield(headmodel, 'cond');
      elseif isfield(headmodel, 'cond') && isfield(headmodel, 'c') && isequal(headmodel.cond, headmodel.c)
        headmodel = rmfield(headmodel, 'cond');
      elseif isfield(headmodel, 'cond') && isfield(headmodel, 'c') && ~isequal(headmodel.cond, headmodel.c)
        ft_error('inconsistent specification of conductive properties for %s model', headmodel.type);
      end
    end

    % ensure that the geometrical units are specified
    headmodel = ft_determine_units(headmodel);

  otherwise
    ft_error('converting to version "%s" is not supported', version);
end

author: spm


text:
%WAVELET  1D Wavelet transform with optional singificance testing
%
%   [WAVE,PERIOD,SCALE,COI,DJ, PARAMOUT, K] = contwt(Y,DT,PAD,DJ,S0,J1,MOTHER,PARAM)
%
%   Computes the continuous wavelet transform of the vector Y (length N),
%   with sampling rate DT.
%
%   By default, the Morlet wavelet (k0=6) is used.
%   The wavelet basis is normalized to have total energy=1 at all scales.
%
%
% INPUTS:
%
%    Y = the time series of length N.
%    DT = amount of time between each Y value, i.e. the sampling time.
%
% OUTPUTS:
%
%    WAVE is the WAVELET transform of Y. This is a complex array
%    of dimensions (N,J1+1). FLOAT(WAVE) gives the WAVELET amplitude,
%    ATAN(IMAGINARY(WAVE),FLOAT(WAVE) gives the WAVELET phase.
%    The WAVELET power spectrum is ABS(WAVE)^2.
%    Its units are sigma^2 (the time series variance).
%
%
% OPTIONAL INPUTS:
%
% *** Note *** setting any of the following to -1 will cause the default
%               value to be used.
%
%    PAD = if set to 1 (default is 0), pad time series with enough zeroes to get
%         N up to the next higher power of 2. This prevents wraparound
%         from the end of the time series to the beginning, and also
%         speeds up the FFT's used to do the wavelet transform.
%         This will not eliminate all edge effects (see COI below).
%
%    DJ = the spacing between discrete scales. Default is 0.25.
%         A smaller # will give better scale resolution, but be slower to plot.
%
%    S0 = the smallest scale of the wavelet.  Default is 2*DT.
%
%    J1 = the # of scales minus one. Scales range from S0 up to S0*2^(J1*DJ),
%        to give a total of (J1+1) scales. Default is J1 = (LOG2(N DT/S0))/DJ.
%
%    MOTHER = the mother wavelet function.
%             The choices are 'MORLET', 'PAUL', or 'DOG'
%
%    PARAM = the mother wavelet parameter.
%            For 'MORLET' this is k0 (wavenumber), default is 6.
%            For 'PAUL' this is m (order), default is 4.
%            For 'DOG' this is m (m-th derivative), default is 2.
%
%
% OPTIONAL OUTPUTS:
%
%    PERIOD = the vector of "Fourier" periods (in time units) that corresponds
%           to the SCALEs.
%
%    SCALE = the vector of scale indices, given by S0*2^(j*DJ), j=0...J1
%            where J1+1 is the total # of scales.
%
%    COI = if specified, then return the Cone-of-Influence, which is a vector
%        of N points that contains the maximum period of useful information
%        at that particular time.
%        Periods greater than this are subject to edge effects.
%        This can be used to plot COI lines on a contour plot by doing:
%
%              contour(time,log(period),log(power))
%              plot(time,log(coi),'k')
%    DJ = actual scale spacing used
%
%    PARAMOUT = actual parameter value used for mother wavelet
%
%   K       =  frequencies at which the transform was computed
%----------------------------------------------------------------------------
%   Copyright (C) 1995-2004, Christopher Torrence and Gilbert P. Compo
%
%   This software may be used, copied, or redistributed as long as it is not
%   sold and this copyright notice is reproduced on each copy made. This
%   routine is provided as is without any express or implied warranties
%   whatsoever.
%
% Notice: Please acknowledge the use of the above software in any publications:
%    ``Wavelet software was provided by C. Torrence and G. Compo,
%      and is available at URL: http://paos.colorado.edu/research/wavelets/''.
%
% Reference: Torrence, C. and G. P. Compo, 1998: A Practical Guide to
%            Wavelet Analysis. <I>Bull. Amer. Meteor. Soc.</I>, 79, 61-78.
%
% Please send a copy of such publications to either C. Torrence or G. Compo:
%  Dr. Christopher Torrence               Dr. Gilbert P. Compo
%  Research Systems, Inc.                 Climate Diagnostics Center
%  4990 Pearl East Circle                 325 Broadway R/CDC1
%  Boulder, CO 80301, USA                 Boulder, CO 80305-3328, USA
%  E-mail: chris[AT]rsinc[DOT]com         E-mail: compo[AT]colorado[DOT]edu
%----------------------------------------------------------------------------
function [wave,period,scale,coi, dj, paramout, k] = ...
	contwt(Y,dt,pad,dj,s0,J1,mother,param);

if (nargin < 8), param = -1;, end
if (nargin < 7), mother = -1;, end
if (nargin < 6), J1 = -1;, end
if (nargin < 5), s0 = -1;, end
if (nargin < 4), dj = -1;, end
if (nargin < 3), pad = 0;, end
if (nargin < 2)
	error('Must input a vector Y and sampling time DT')
end

n1 = length(Y);

if (s0 == -1), s0=2*dt;, end
if (dj == -1), dj = 1./4.;, end
if (J1 == -1), J1=fix((log(n1*dt/s0)/log(2))/dj);, end
if (mother == -1), mother = 'MORLET';, end

%....construct time series to analyze, pad if necessary
x(1:n1) = Y - mean(Y);
if (pad == 1)
	base2 = fix(log(n1)/log(2) + 0.4999);   % power of 2 nearest to N
	x = [x,zeros(1,2^(base2+1)-n1)];
end
n = length(x);

%....construct wavenumber array used in transform [Eqn(5)]
k = [1:fix(n/2)];
k = k.*((2.*pi)/(n*dt));
k = [0., k, -k(fix((n-1)/2):-1:1)];

%....compute FFT of the (padded) time series
f = fft(x);    % [Eqn(3)]

%....construct SCALE array & empty PERIOD & WAVE arrays
scale = s0*2.^((0:J1)*dj);
period = scale;
wave = zeros(J1+1,n);  % define the wavelet array
wave = wave + i*wave;  % make it complex

% loop through all scales and compute transform
for a1 = 1:J1+1
	[daughter,fourier_factor,coi,dofmin, paramout]=wave_bases(mother,k,scale(a1),param);
	wave(a1,:) = ifft(f.*daughter);  % wavelet transform[Eqn(4)]
end

period = fourier_factor*scale;
coi = coi*dt*[1E-5,1:((n1+1)/2-1),fliplr((1:(n1/2-1))),1E-5];  % COI [Sec.3g]
wave = wave(:,1:n1);  % get rid of padding before returning

return

% end of code

author: Sable

repo: mcbench-benchmarks


text:
function newCoords = UpsampleCoords(coords, factor)

% newCoords = UpsampleCoords(coords, factor);
%
% Returns a coordinate array that oversamples the input coords by the
% specified factor.
%
% Ress, 2/04

% Determine a few index parameters
factor = ceil(factor);
i1 = factor - 1;
i0 = -i1;
stride = 1 / factor;
dims = size(coords);
nCoords = dims(2);

% Initialize the upsampled coordinates:
newCoords = zeros(3, nCoords*(2*factor-1)^3);
iC = 1;

for iz=i0:i1
  offset(3) = iz * stride;
  for iy=i0:i1
    offset(2) = iy * stride;
    for ix=i0:i1
      offset(1) = ix * stride;
      offsetCoords = coords;
      for ii=1:3, offsetCoords(ii, :) = offsetCoords(ii, :) + offset(ii); end
      newCoords(:, iC:(iC+nCoords-1)) = offsetCoords;
      iC = iC + nCoords;
    end
  end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function preds = predict_box(box, modelPath, saveConfmaps)
%PREDICT_BOX Evaluates model predictions on a stack of frames. Wrapper for predict_box.py.
% Usage:
%   preds = predict_box(box, modelPath)
%   preds = predict_box(box, modelPath, saveConfmaps)
%
% Args:
%   box: 4-D array or path to HDF5 file with '/box' dataset
%   modelPath: path to model weights
%   saveConfmaps: if true, returns full confidence maps in addition to
%                 peaks (default: false). Very slow and memory intensive!
%
% Returns:
%   preds: struct containing results from model prediction
%        .positions_pred: 3-D array of (parts x [X Y] x frames) indicating
%                         peak positions for each confidence map in image coordinates
%        .conf_pred: 2-D array of (parts x frames) with the confidence map
%                    value at the peak pixel for detecting bad predictions
%        .confmaps: 4-D array of confidence maps returned if saveConfmaps
%                   is set true

if nargin < 3 || isempty(saveConfmaps); saveConfmaps = false; end

% Process args
delete_box = false;
is_singleton = false;
if ischar(box)
    boxPath = box;
else
    boxPath = [tempname '.h5'];

    if numel(size(box)) < 4
        box = repmat(box,[1 1 1 2]);
        is_singleton = true;
    end

    h5save(boxPath, box)
    delete_box = true;
end

% Generate temporary output filename
outPath = [tempname '.h5'];

% Build command line args
cmd = {
    'python'
    ['"' ff(funpath(true), 'predict_box.py') '"']
    ['"' boxPath '"']
    ['"' modelPath '"']
    ['"' outPath '"']
    };
if saveConfmaps
    cmd{end+1} = '--save-confmaps';
end
disp(strjoin(cmd))

% Predict
try
    exit_code = system(strjoin(cmd));
catch ME
    if delete_box && exists(boxPath); delete(boxPath); end
    rethrow(ME)
end
if delete_box && exists(boxPath); delete(boxPath); end

% Read data back in
try
    preds = h5readgroup(outPath);
catch ME
    if exists(outPath); delete(outPath); end
    rethrow(ME)
end
if exists(outPath); delete(outPath); end

% Adjust for 0-based indexing
preds.positions_pred = single(preds.positions_pred) + 1;

% Rescale confidence maps to correct range prior to quantization
if saveConfmaps && isfield(preds.Attributes, 'confmaps')
    c_min = preds.Attributes.confmaps.range_min;
    c_max = preds.Attributes.confmaps.range_max;
    preds.confmaps = rescale(single(preds.confmaps) / 255, c_min, c_max);
end

% Adjust for singleton input
if is_singleton
    preds.positions_pred = preds.positions_pred(:,:,1);
    preds.conf_pred = preds.conf_pred(:,1);
    if isfield(preds,'confmaps')
        preds.confmaps = preds.confmaps(:,:,:,1);
    end
end

end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
function [ t_stop, y_stop ] = p17_stop ( neqn )

%*****************************************************************************80
%
%% P17_STOP returns the stopping point for problem p17.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_STOP, Y_STOP(NEQN), the final data.
%
  y_stop = zeros ( neqn, 1 );

  t_stop = 20.0;
  y_stop = [ ...
   -1.777027357140412e-01; ...
    9.467784719905892e-01; ...
   -1.030294163192969; ...
    1.211074890053952e-01 ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [K, R, t] = decomposeP(P)
%VGG_KR_FROM_P Extract K, R from camera matrix.
%
%    [K,R,t] = VGG_KR_FROM_P(P [,noscale]) finds K, R, t such that P = K*R*[eye(3) -t].
%    It is det(R)==1.
%    K is scaled so that K(3,3)==1 and K(1,1)>0. Optional parameter noscale prevents this.
%
%    Works also generally for any P of size N-by-(N+1).
%    Works also for P of size N-by-N, then t is not computed.

% Author: Andrew Fitzgibbon <awf@robots.ox.ac.uk>
% Modified by werner.
% Date: 15 May 98

N = size(P,1);
H = P(:,1:N);

[K,R] = rq(H);

K = K*sign(K(1));
if prod(sign(diag(K))) < 0
  D = diag(-1*sign(diag(K)));
else
  D = diag(sign(diag(K)));
end
K = K*D;
R = D*R;

K = K/K(end);

% $$$ if nargin < 2
% $$$   K = K / K(N,N);
% $$$   if K(1,1) < 0
% $$$     D = diag([-1 -1 ones(1,N-2)]);
% $$$     K = K * D;
% $$$     R = D * R;
% $$$
% $$$   %  test = K*R;
% $$$   %  vgg_assert0(test/test(1,1) - H/H(1,1), 1e-07)
% $$$   end
% $$$ end

if nargout > 2
  t = -P(:,1:N)\P(:,end);
end

return

author: jianxiongxiao

repo: ProfXkit


text:
                    function rot
% ***************************************************************
%                   P r o g r a m   ROT
% ***************************************************************
%
%  PURPOSE:
%     Resolve numerically differential equation of rotational motion
%     of a body
%                      Jz*d2f/dt2 = Mz(t,f,w)
%     and plots graphics of coordinate, velocity and phase plane.
%     If possible, the program could solve the problem analytically.
%
%  INPUT DATA:
%     Jz   - moment of inertia of the body ;
%     Mz   - rotational moment  Mz = Mz(t,f,w);
%     f0   - initial value of the coordinate ;
%     w0   - initial value of the angular velocity ;
%     Tend - upper bound of the integration ;
%     eps  - precision of the integration ;
%     np   - number of parameters .
%     P{1}, P{2}, ..., P{np} - names of the parameters (array of cells);
%
%  NOTES:
%   1. The coordinate is designed by the symbol 'f' and velocity by 'w';
%   2. The physical names of the parameters are assigned to the
%      cells of the array P like this: P{1}='Jz', P{2}='c',...;
%   3. For analytical solution the values of Tend, eps, np and P are not
%      needed.
%   4. The parameters Jz, f0 and w0 have to be entered only as
%      strings, even though they represent numbers!
%   5. All the data can be input from file or in interactive mode.
%
%  EXAMPLE of DATA FILE:
%   % Data File for problem ...
%     Jz   = 'Jz'; ( or Jz = '0.15';)
%     Mz   = '-k*w - c*f'; % k, c - parameters
%     f0   = 'f0'; ( or f0 = '0.33';)
%     w0   = 'w0'; ( or w0 = '7';)
%     Tend = 20;
%     eps  = 1.e-8;
%     np   = 3;
%     P{1} = 'Jz';
%     P{2} = 'k';
%     P{3} = 'c';

% ---------------------------------------------------------
%               DATA INPUT OF THE PROBLEM
% =========================================================

 clear
 disp(' ');
 disp(' How will you input the data ?    ');
 disp('     1. From a data file;         ');
 disp('     2. In interactive mode.      ');
 ans = input(' Number of Your choice : '  );
 flag = 0;
 if ans == 1
    while 1
       disp(' ');
       indat = input(' Input the name of the data file :', 's');
       if exist([cd,'\',indat]) % Search only in current directory
          eval(indat);
          flag = 1; break  % Successful
        else               % Unsuccessful
          disp(' ');
          disp([' File ',indat,' not exist!'])
          disp(' You have to:')
          disp(' 1. Enter another DATA file name, or')
          disp(' 2. Input the DATA interactively !')
          ans2 = input(' Your choice, please: ');
          if ans2 == 2, break , end
       end
    end
 end
if flag == 0
    %Input of the data in-line mode
    disp(' ');
    Jz   = input(' Moment of inertia of the body Jz : ','s'      );
    Mz   = input(' Expression of the rotational moment Mz : ','s');
    f0   = input(' Initial coordinate f0 : '                     );
    w0   = input(' Initial velocity w0 : '                       );
    Tend = input(' Upper bound of the integration Tend : '       );
    eps  = input(' Precision of the calculations eps : '         );
    np   = input(' Number of parameters  np : '                  );
    % Asigning names of the parameters
    if np > 0
       disp(' ');
       disp(' Enter names of the parameters:')
       for i = 1:np
           ii = num2str(i);
           P{i} = input([' Name of parameter P',ii,': '],'s');
       end
    end
 end

% ---------------------------------------------------------
%            Differential Equation of Motion
% =========================================================

                 syms f w D2f
                 Mz = subs(Mz, 'w', 'Df');
                 deq = Jz*D2f - Mz;

% ---------------------------------------------------------
%                  ANALYTICAL SOLUTION
% =========================================================

 disp(' ');
 ans = input(' Would you like analytical solution? (Y/N): ','s');
 if ans =='Y' | ans == 'y'
    if ~isstr(f0), f0 = num2str(f0); end % Repairing user
    if ~isstr(w0), w0 = num2str(w0); end % input errors!
    syms f
    inicond = ['f(0)=',f0,',Df(0)=',w0];
    f = dsolve(char(deq), inicond, 't');
    if ~isempty(f)
        disp(' ');
        disp('   Low of Motion   ');
        disp(' ***************** ');
        disp(' ');
        disp('f = '); pretty(f)
        disp(' ');
        fname = input(' Name of file to write solution: ','s');
        save(fname, 'f');
    end
       disp(' ');
       ans = input(' Would you like numerical solution? (Y/N): ','s');
       if ans == 'N' | ans == 'n', return, end
       f = 'f'; % Clear contents of f
 end

% ---------------------------------------------------------
%                  NUMERICAL SOLUTION
% =========================================================

% Input the name of the file-function
disp(' ');
fname = input(' Name of the File-function to be generated: ','s');
flag1 = 'Y';
if exist([cd,'\',fname]) % Search only in current directory!
    disp(' ');
    disp([' A file-function with name ',fname,' already exist !'])
    flag1 = input(' Overwrite it ? (Y/N): ', 's');
end

% ---------------------------------------------------------
%              GENERATING THE FILE-FUNCTION
% ---------------------------------------------------------

if ( flag1 == 'Y' | flag1 == 'y' )
   Mz = subs(Mz,{'f','Df'},{'y(1)','y(2)'});
% Opening the file for writing file-function
   [Fid,mes] = fopen([fname,'.m'],'wt');
% Generating the string with physical parameters: Jz, c ...
   strpar = '';
   for j = 1:np
      strpar = [strpar,',',P{j}];
   end
   titl = input(' Denomination of the Problem: ','s');
%       Writing the headline of the File-function
   fprintf(Fid,['function yt = ',fname,'(t,y',strpar,')\n']);
   fprintf(Fid,['%% ',titl]);
%       Writing the first derivatives
   fprintf(Fid,'\n%% The first derivatives\n');
   fprintf(Fid, '  yt(1) = y(2); \n');
   fprintf(Fid,['  yt(2) = ',char(Mz/Jz),'; \n']);
   fprintf(Fid,'  yt = yt'';\n');
   fprintf(Fid,['%% *** End of File-function ',fname,' ***']);
   fclose(Fid);
   edit(fname)
end

% ---------------------------------------------------------
%        INTEGRATION AND VISUALIZATION OF THE REZULTS
% ---------------------------------------------------------

flag2 = 0;
% Initial entering values of the parameters and generating
% the string with parameters 'P{1}, P{2}, ..., P{np}' to be
% passed to the File-function as actual arguments
if np > 0
   PP = P; % Saving the physical names of the parameters in PP
   parameters = ' ';
   disp(' ');
   disp(' Input the numerical values of the parameters: ')
   for i = 1:np
       i = num2str(i);
       eval(['P{',i,'}=input([''   '',P{',i,'},'' = '']);']);
       parameters = [parameters,',P{',i,'}'];
   end
else
    parameters = [];
end
% Check-up type of f0 and w0 and correct
% it if needed
if ischar(f0)
    f0 = str2num(f0);
    if isempty(f0), f0 = input(' f0 = '); end
end
if ischar(w0)
    w0 = str2num(w0);
    if isempty(w0), w0 = input(' w0 = '); end
end
while 1
    if flag2 == 1
        disp(' ');
        eps  = input(' Precision of the computations eps: ');
        Tend = input(' Upper bound of the integration Tend: ');
        f0   = input(' Initial coordinate f0: ');
        w0   = input(' Initial velocity w0: ');
        if np > 0
          P = PP; % Restoring the names of the parameters !
          disp(' ');
          disp(' Input the numerical values of the parameters: ')
          for i = 1:np
              i = num2str(i);
              eval(['P{',i,'}=input([''   '',P{',i,'},'' = '']);']);
          end
        end
    end
    y0 = [f0 w0]; % initial conditions
    options = odeset('AbsTol',eps,'RelTol',100*eps);
    % Choosing of the Solver
    disp('                                        ');
    disp('      Choose the proper Solver:         ');
    disp('  -------------------------------       ');
    disp(' A. Non stiff differential equations    ');
    disp('   1. ode45   - middle precision;       ');
    disp('   2. ode23   - low precision;          ');
    disp('   3. ode113  - from low to upper.      ');
    disp('                                        ');
    disp(' B. Stiff differential equations        ');
    disp('   1. ode15s  - from low to upper;      ');
    disp('   2. ode23s  - low precision;          ');
    disp('   3. ode23t  - middle precision;       ');
    disp('   4. ode23tb - low precision.          ');
    disp('                                        ');
    solver = input(' The name of the Solver: ','s');

    % Integration of the Differential Equations

    eval(['[t,y] = feval(solver,eval([''@'',fname]),',...
                  '[0 Tend],y0,options',parameters,');']);
    % Plotting graphs

    tmin  = min(t);
    tmax  = max(t);
    y1min = min(y(:,1));
    y1max = max(y(:,1));
    y2min = min(y(:,2));
    y2max = max(y(:,2));
    dy1   = y1max - y1min;
    dy2   = y2max - y2min;
    xmin  = y1min - 0.1*dy1;
    xmax  = y1max + 0.1*dy1;
    ymin  = y2min - 0.1*dy2;
    ymax  = y2max + 0.1*dy2;
    % Coordinate f = f(t)
    figure % 1
    comet(t,y(:,1))
    plot(t,y(:,1),[tmin tmax],[0 0],'k'), grid on
    axis([tmin, tmax, xmin, xmax]);
    set(gca,'FontName','Arial Cyr','FontSize',12);
    title('Low of motion {\phi} = {\phi}({\itt})')
    xlabel('{\itt}'); ylabel('{\phi}'); pause
    % Angular velocity w = w(t)
    figure % 2
    comet(t,y(:,2))
    plot(t,y(:,2),[tmin tmax],[0 0],'k'), grid on
    axis([tmin, tmax, ymin, ymax]);
    set(gca,'FontName','Arial','FontSize',12);
    title('Angular velocity {\omega} = {\omega}({\itt})')
    xlabel('{\itt}'); ylabel('{\omega}'); pause
    % Coordinate and Velocity
    figure % 3
    subplot(2,1,1), plot(t,y(:,1),[tmin tmax],[0 0],'k')
    grid on, axis([tmin, tmax, xmin, xmax]);
    set(gca,'FontName','Arial','FontSize',12);
    title('Low of motion {\phi} = {\phi}(t)')
    subplot(2,1,2), plot(t,y(:,2),[tmin tmax],[0 0],'k')
    grid on, axis([tmin, tmax, ymin, ymax]);
    set(gca,'FontName','Arial','FontSize',12);
    title('Angular velocity {\omega} = {\omega}(t)')
    pause
    % Phase Plane
    figure % 4
    subplot(1,1,1)
    comet(y(:,1),y(:,2))
    plot(y(:,1),y(:,2), [xmin,xmax],[0 0],'k',...
                  [0 0],[ymin,ymax],'k'), grid on
    axis([xmin, xmax, ymin, ymax]);
    set(gca,'FontName','Arial Cyr','FontSize',12);
    title(' Phase Plane {\omega} = {\omega}({\phi})')
    xlabel('{\phi}'), ylabel('{\omega}'), pause
    flag2 = 1;
    % close all
    disp(' ');
    ans = input(' Would you like to continue? (Y/N): ','s');
    if ans == 'n' | ans == 'N', break, end
end

%  ***************** End of Program ROT ********************


text:
function varargout = gui_stat_analyticStats(varargin)
%
% GUI_STAT_ANALYTICSTATS M-file for gui_stat_analyticStats.fig
%      GUI_STAT_ANALYTICSTATS, by itself, creates a new GUI_STAT_ANALYTICSTATS or raises the existing
%      singleton*.
%
%      H = GUI_STAT_ANALYTICSTATS returns the handle to a new GUI_STAT_ANALYTICSTATS or the handle to
%      the existing singleton*.
%
%      GUI_STAT_ANALYTICSTATS('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in GUI_STAT_ANALYTICSTATS.M with the given input arguments.
%
%      GUI_STAT_ANALYTICSTATS('Property','Value',...) creates a new GUI_STAT_ANALYTICSTATS or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before gui_stat_analyticStats_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to gui_stat_analyticStats_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help gui_stat_analyticStats

% Last Modified by GUIDE v2.5 28-Jun-2012 15:12:51

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @gui_stat_analyticStats_OpeningFcn, ...
                   'gui_OutputFcn',  @gui_stat_analyticStats_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before gui_stat_analyticStats is made visible.
function gui_stat_analyticStats_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to gui_stat_analyticStats (see VARARGIN)

set(hObject,'name','Analytic Statistics');

handles.output = hObject;

% set default termination behavior
handles.ExitButtonClicked = 'Cancel';

% extract some data from command-line input
if isempty(varargin)
    error('You must pass ALLEEG to gui_analyticStats');
end

% Extract input parameters/data and store
handles.ud.ALLEEG  = varargin{1};
varargin(1) = [];

% set default EEGLAB background and text colors
%-----------------------------------------------
try, icadefs;
catch,
	GUIBACKCOLOR        =   [.8 .8 .8];
	GUIPOPBUTTONCOLOR   =   [.8 .8 .8];
	GUITEXTCOLOR        =   [0 0 0];
end;

allhandlers = hObject;

hh = findobj(allhandlers,'style', 'text');
%set(hh, 'BackgroundColor', get(hObject, 'color'), 'horizontalalignment', 'left');
set(hh, 'Backgroundcolor', GUIBACKCOLOR);
set(hh, 'foregroundcolor', GUITEXTCOLOR);
set(hObject, 'color',GUIBACKCOLOR );
% set(hh, 'horizontalalignment', g.horizontalalignment);

hh = findobj(allhandlers, 'style', 'edit');
set(hh, 'BackgroundColor', [1 1 1]); %, 'horizontalalignment', 'right');

hh =findobj(allhandlers, 'parent', hObject, 'style', 'pushbutton');
if ~strcmpi(computer, 'MAC') && ~strcmpi(computer, 'MACI') % this puts the wrong background on macs
    set(hh, 'backgroundcolor', GUIPOPBUTTONCOLOR);
    set(hh, 'foregroundcolor', GUITEXTCOLOR);
end;
hh =findobj(allhandlers, 'parent', hObject, 'style', 'popupmenu');
set(hh, 'backgroundcolor', GUIPOPBUTTONCOLOR);
set(hh, 'foregroundcolor', GUITEXTCOLOR);
hh =findobj(allhandlers, 'parent', hObject, 'style', 'checkbox');
set(hh, 'backgroundcolor', GUIBACKCOLOR);
set(hh, 'foregroundcolor', GUITEXTCOLOR);
hh =findobj(allhandlers, 'parent', hObject, 'style', 'listbox');
set(hh, 'backgroundcolor', GUIPOPBUTTONCOLOR);
set(hh, 'foregroundcolor', GUITEXTCOLOR);
hh =findobj(allhandlers, 'parent', hObject, 'style', 'radio');
set(hh, 'foregroundcolor', GUITEXTCOLOR);
set(hh, 'backgroundcolor', GUIPOPBUTTONCOLOR);
set(hObject, 'visible', 'on');

set(handles.pnlPropertyGrid,'backgroundcolor', GUIBACKCOLOR);
set(handles.pnlPropertyGrid,'foregroundcolor', GUITEXTCOLOR);
%-----------------------------------------------

drawnow

% render the PropertyGrid in the correct panel
handles.PropertyGridHandle = arg_guipanel( ...
                 handles.pnlPropertyGrid, ...
                'Function',@stat_analyticStats, ...
                'params',{'ALLEEG',handles.ud.ALLEEG, varargin{:}});

% Update handles structure
guidata(hObject, handles);

% Wait for user to click OK, Cancel or close figure
uiwait(handles.gui_analyticStats);

% --- Outputs from this function are returned to the command line.
function varargout = gui_stat_analyticStats_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isempty(handles)
    % user closed the figure
    varargout = {[] hObject};
elseif strcmpi(handles.ExitButtonClicked,'OK')
    % user clicked OK
    % get PropertySpecification
    varargout = {handles.PropertyGridHandle handles.output};
else
    % user clicked cancel
    varargout = {[] handles.output};
end

try, close(hObject);
catch; end

function cmdCancel_Callback(hObject, eventdata, handles)

handles.ExitButtonClicked = 'Cancel';
guidata(hObject,handles);
uiresume(handles.gui_analyticStats)

function cmdOK_Callback(hObject, eventdata, handles)

handles.ExitButtonClicked ='OK';
guidata(hObject,handles);
uiresume(handles.gui_analyticStats);

function cmdHelp_Callback(hObject, eventdata, handles)
doc('stat_analyticStats');

function gui_analyticStats_ResizeFcn(hObject, eventdata, handles)

repo: REST


text:
function sz = mysize(M)
% MYSIZE Like the built-in size, except it returns n if M is a vector of length n, and 1 if M is a scalar.
% sz = mysize(M)
%
% The behavior is best explained by examples
% - M = rand(1,1),   mysize(M) = 1,      size(M) = [1 1]
% - M = rand(2,1),   mysize(M) = 2,      size(M) = [2 1]
% - M = rand(1,2),   mysize(M) = 2,      size(M) = [1 2]
% - M = rand(2,2,1), mysize(M) = [2 2],  size(M) = [2 2]
% - M = rand(1,2,1), mysize(M) = 2,      size(M) = [1 2]

if isvectorBNT(M)
  sz = length(M);
else
  sz = size(M);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function au_qmr_ilu_s_i(p)
%
%  Generates the data used in 'au_qmr_ilu_s'. Data are stored in global
%  variables.
%
%  NOTE that 'au_qmr_ilu_m_i' must be called before calling this routine.
%
%  Calling sequence:
%
%    au_qmr_ilu_s_i(p)
%
%  Input:
%
%    p         vector containing the ADI shift parameters p(i).
%
%  Remarks:
%
%    This routine has access to the matrix A and some other data, which
%    must be provided by the routine 'au_qmr_ilu_m_i',
%
%    The real parts of the entries of p must be negative.
%
%
%  LYAPACK 1.0 (Thilo Penzl, August 1999)

if nargin~=1
  error('Wrong number of input arguments.');
end

if any(real(p)>=0)
  error('Real parts of entries of p must be negative!');
end

l = length(p);

global LP_A LP_P LP_MC LP_TOL_ILU LP_INFO_QMR

if ~length(LP_A) | ~length(LP_MC) | ~length(LP_TOL_ILU) | ~length(LP_INFO_QMR)
  error('This routine needs global data which must be generated by calling ''au_qmr_ilu_m_i'' first.');
end

n = size(LP_A,1);

LP_P = p;

if LP_MC=='C'

  for i = 1:l

    eval(lp_e( 'global LP_L',i,' LP_U',i ));

    eval(lp_e( '[LP_L',i,',LP_U',i,'] = luinc(LP_A+p(i)*speye(n),LP_TOL_ILU);' ));

%
%  This shows the amount of fill-in created by the ILU preconditioner:
%

    if LP_INFO_QMR >= 5
      nnA = nnz(LP_A);
      eval(lp_e( 'nnLU = nnz(LP_L',i,')+nnz(LP_U',i,');'  ));
      disp(sprintf('ILU-Ratio (nnz(L)+nnz(U))/nnz(A) = %6g',nnLU/nnA));
    end
  end

end


text:
function updateSpatial_endoscope(obj, Y, num, method, smin)
%% udpate spatial components

%% inputs:
%   Y: d X T matrix, data
%   num: scalar. If method=='hals', then num is the number of iterations to
%       update A; If method=='nnls', then num is the maximum number of neurons
%       overlapping at one pixel
%   method: method for updating the spatial components {'hals', 'nnls'}.
%       default: 'nnls'

%% Author: Pengcheng Zhou, Carnegie Mellon University.

%% input parameters number of iterations
if ~exist('method', 'var')||isempty(method)
    method = 'nnls';
end
if ~exist('num', 'var')||isempty(num)
    if strcmpi(method, 'nnls')
        num=5;
    else
        num = 10;
    end
end
if ~exist('IND_thresh', 'var')||isempty(IND_thresh)
    IND_thresh = [];
end
%% determine the search locations
search_method = obj.options.search_method;
params = obj.options;
if strcmpi(search_method, 'dilate')
    obj.options.se = [];
end
IND = logical(determine_search_location(obj.A, search_method, params));

%% update spatial components
if strcmpi(method, 'hals')
    obj.A = HALS_spatial(Y, obj.A, obj.C, IND, num);
elseif strcmpi(method, 'nnls_thresh')&&(~isempty(IND_thresh))
    try
        sn = obj.P.sn;
    catch
        sn = get_noise_fft(Y);
        obj.P.sn = sn;
    end

    obj.A = nnls_spatial_thresh(Y, obj.A, obj.C, IND, num, smin, sn);
else
    obj.A = nnls_spatial(Y, obj.A, obj.C, IND, num);
end

%% thresholding the minimum number of neurons
obj.delete(sum(obj.A, 1)<=obj.options.min_pixel);

%% post-process
% obj.post_process_spatial();
% if strcmpi(method, 'nnls')
%     IND = bsxfun(@gt, obj.A, max(obj.A, [], 1)/100);
%     obj.A = nnls_spatial(Y, obj.A, obj.C, IND, num);
% end

end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function [] = savefig(image_name,font_size,formats)
%   SAVEFIG   save current figure in image file
%       [] = SAVEFIG(IMAGE_NAME,FONT_SIZE)
%
%   Created by Alexandre Gramfort on 2008-11-11.
%   Copyright (c) 2007 Alexandre Gramfort. All rights reserved.

% $Id: savefig.m 4 2009-08-15 21:10:35Z gramfort $
% $LastChangedBy: gramfort $
% $LastChangedDate: 2009-08-15 17:10:35 -0400 (Sam, 15 aoû 2009) $
% $Revision: 4 $

% hack
return

if nargin < 3
    formats = {'png'};
end

if nargin < 2
    font_size = 22;
end

if ~exist('image_name')
    image_name = get(gcf,'Name');
end

all_axis = get(gcf,'children');

name_eps = [image_name,'.eps'];
name_jpg = [image_name,'.jpg'];
name_png = [image_name,'.png'];
name_pdf = [image_name,'.pdf'];

for ii=1:length(all_axis)
    try
        title(all_axis(ii),' ')
        set(all_axis(ii),'FontSize',font_size);
        set(get(all_axis(ii),'xlabel'),'FontSize',font_size);
        set(get(all_axis(ii),'ylabel'),'FontSize',font_size);
        set(get(all_axis(ii),'zlabel'),'FontSize',font_size);
    end
end

if any(strcmp('png',formats))
    saveas(gcf,name_png,'png');
    %%% try to trim png file
    disp('please wait, trim png file....');
    command = sprintf(['convert -trim %s %s'],name_png,name_png);
    [STATUSVAR,STATUSMESSAGE] = unix(command);

    if STATUSVAR ~= 0
       error( STATUSMESSAGE );
    else
       disp('trim completed');
    end
end

if any(strcmp('pdf',formats))
    print_pdf(name_pdf,gcf);
end

if any(strcmp('eps',formats))
    saveas(gcf,name_eps,'epsc');
end

if any(strcmp('jpg',formats))
    saveas(gcf,name_jpg,'jpg');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Prepare a network for 'test' mode
% - disable `precious`
% - optionally remove `dropout`

function net= netPrepareForTest(net, removeDropout, removeExtras)
    if nargin<2, removeDropout= false; end
    if nargin<3, removeExtras= false; end

    for iLayer= 1:length(net.layers)
        net.layers{iLayer}.precious= false;

        if removeExtras
            if isFieldOrProp(net.layers{iLayer}, 'momentum')
                net.layers{iLayer}.momentum= [];
            end
        end
    end

    if removeDropout
        net.layers( ismember( relja_layerTypes(net), {'dropout'} ) )= [];
    end

end

function is= isFieldOrProp(l, propName)
    is= isprop(l, propName) || isfield(l, propName);
end

author: Relja

repo: netvlad

save_path: github-repos/MATLAB/Relja-netvlad


text:
% MAIN.m -- Robot Arm Trajectory
%
% Given an initial and final state, as well as joint angle and rate
% constraints, find the minimum-jerk trajectory that solves the boundary
% value problem.
%
%

n = 15;  %Order of interpolating polynomial;

% Problem specifications:
tStart = 0;
tFinal = 1;
qLow = -pi;
qUpp = pi;
dqMax = 1.5;   %Joint speed limit

% Boundary values
qStart = 0;
qFinal = 1;
dqStart = 0;
dqFinal = -1.5;

% Derivatives
D = chebyshevDifferentiationMatrix(n,[tStart,tFinal]);   %Rates
DD = D*D;    %Accelerations
DDD = DD*D;   %Jerk

% Chebyshev grid and quadrature weights:
[t,w] = chebyshevPoints(n,[tStart,tFinal]);

%%%% Quadratic Program:
%   min 0.5*x'*H*x + f'*x   subject to:  A*x <= b
%    x                                   Aeq*x == beq
%                                        LB <= x <=UB

%%%% Boundary Values:
Aeq = zeros(4,n); beq = zeros(4,1);
Aeq(1,1) = 1; beq(1) = qStart;
Aeq(2,end) = 1; beq(2) = qFinal;
Aeq(3:4,:) = D([1,end],:); beq(3:4) = [dqStart; dqFinal];

%%%% Dynamics:
% u == DDD*x;   %Equality Constraint - expressed directly in cost function

%%%% Joint Rate Limits:
% -dqMax <= D*x <= dqMax
A = [D;-D];
b = [dqMax*ones(n,1); dqMax*ones(n,1)];

%%%% Joint Angle Limits:
% qLow <= x <= qUpp
LB = qLow*ones(n,1);
UB = qUpp*ones(n,1);

%%%% Objective Function:
% min  w*(u^2) = x'*(DDD'*W*DDD);  W = diag(w);
W = diag(w);
H = 0.5*(DDD'*W*DDD);

%%%% Initialize:    (Linear trajectory)
xGuess = interp1([tStart;tFinal],[qStart;qFinal],t');

%%%% Options:
options = optimset(...
    'Display','iter',...  % {'iter','final'}
    'Algorithm','interior-point-convex');

%%%% Build Problem:
problem.H = H;
problem.f = zeros(n,1);
problem.Aineq = A;
problem.bineq = b;
problem.Aeq = Aeq;
problem.beq = beq;
problem.lb = LB;
problem.ub = UB;
problem.x0 = xGuess;
problem.options = options;
problem.solver = 'quadprog';

%%%% Solve!
[x, fVal, exitFlag, output] = quadprog(problem);

% Plot the solution:
tt = linspace(t(1),t(end),100);
[xx, dxx] = chebyshevInterpolate(x',tt,[t(1),t(end)]);
figure(2); clf;
subplot(2,1,1); hold on;
plot(tt,xx);
plot(t,x,'ko');
subplot(2,1,2); hold on;
plot(tt,dxx);
plot(t,D*x,'ko');

author: MatthewPeterKelly

repo: dscTutorials


text:
function L = spm_mesh_get_lm(M,T)
% Identification of local maxima on a textured surface mesh
% FORMAT L = spm_mesh_get_lm(M,T)
% M        - a [nx3] faces array or a patch structure or a [nxn] adjacency
%            matrix
% T        - a [nx1] texture vector
%
% L        - indices of vertices that are local maxima
%__________________________________________________________________________
% Copyright (C) 2010-2016 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_mesh_get_lm.m 6867 2016-09-12 15:04:44Z guillaume $

%-Get adjacency matrix
%--------------------------------------------------------------------------
if ~isnumeric(M) || size(M,1) ~= size(M,2)
    A = spm_mesh_adjacency(M);
else
    A = M;
end

%-Get neighbours, restricted to vertices that actually have data (~= NaN)
%--------------------------------------------------------------------------
out      = isnan(T(:)');
if all(out), L = []; return; end             % empty domain
A(:,out) = 0;
A(out,:) = 0;
N        = spm_mesh_neighbours(A);
if isempty(N), L = find(~out); return; end   % only singletons
%S       = ~any(N(~out,:),2);                % some singletons

%-Identify local maxima
%--------------------------------------------------------------------------
T        = [T; -Inf];
N(N<1)   = numel(T);
L        = all(bsxfun(@gt,T(1:end-1),T(N)),2);
L        = find(L');


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nickabattista@gmail.com
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs, non-invariant beams*)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "Hill+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific muscle model, please let Nick (nickabattista@gmail.com) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: updates the spring attributes!
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function springs_info = update_Springs(dt,current_time,xLag,yLag,springs_info)

%springs_info: col 1: starting spring pt (by lag. discretization)
%              col 2: ending spring pt. (by lag. discretization)
%              col 3: spring stiffness
%              col 4: spring resting lengths

%RL = springs_info(:,4); % resting-length vector

% Contraction Frequency
freq = 1.0;

% CHANGE RESTING LENGTH BTWN SIDES OF JELLYFISH BELL
springs_info(637:end,4) = abs( cos(freq*current_time*pi) );

%NOTE: not 2*pi*ft b/c of abs()

author: nickabattista


text:
classdef PoseNode < handle
    %POSENODE A class for pose node

    properties (Access = public)
        id    % Id of this pose node
        pose  % Pose of this pose node
    end  % properties public

    properties (Dependent = true)
        x    % X coordinate
        y    % Y coordinate
        yaw  % Yaw angle
        rt   % Transformation local to global
    end  % properties dependent

    methods

        function obj = PoseNode(id, pose)
            % Constructor of PoseNode
            obj.id   = id;
            obj.pose = pose(:);
        end

        function plot(obj)
            % Plot all pose nodes position
            x = [obj.x];
            y = [obj.y];
            plot(x, y, 'b');
        end

        function x = get.x(obj)
            x = obj.pose(1);
        end

        function y = get.y(obj)
            y = obj.pose(2);
        end

        function yaw = get.yaw(obj)
            yaw = obj.pose(3);
        end

        function rt = get.rt(obj)
            R = [cos(obj.yaw) -sin(obj.yaw);
                 sin(obj.yaw)  cos(obj.yaw)];
            rt = [R [obj.x; obj.y]; 0 0 1];
        end

    end  % methods public

end  % classdef

repo: graphslam


text:
% SP_HCURL_ERROR: Evaluate the error in H(curl) norm.
%
%   [errhcurl, errl2, errcurl] = sp_hcurl_error (space, msh, u, uex, curluex);
%
% INPUT:
%
%    space:   struct defining the space of discrete functions (see sp_vector/sp_evaluate_col)
%    msh:     struct defining the domain partition and the quadrature rule (see msh_cartesian/msh_evaluate_col)
%    u:       vector of dof weights
%    uex:     function handle to evaluate the exact solution
%    curluex: function handle to evaluate the curl of the exact solution
%
% OUTPUT:
%
%     errhcurl: error in H(curl) norm
%     errl2:    error in L^2 norm
%     errcurl:  error of the curl in L^2 norm
%
% Copyright (C) 2010 Carlo de Falco, Rafael Vazquez
% Copyright (C) 2011, 2015 Rafael Vazquez
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Octave; see the file COPYING.  If not, see
% <http://www.gnu.org/licenses/>.

function [errhcurl, errl2, errcurl] = sp_hcurl_error (sp, msh, u, uex, curluex)

  w = msh.quad_weights(:) .* msh.jacdet(:);

  for idim = 1:msh.rdim
    x{idim} = reshape (msh.geo_map(idim,:,:), msh.nqn*msh.nel, 1);
  end

  curl_valex = feval (curluex, x{:});
  errcurl = 0;

  switch (msh.rdim)
   case {2}
     valnum = zeros (msh.nqn, msh.nel);
     for ish = 1:sp.nsh_max
       valnum = valnum + ...
         reshape (sp.shape_function_curls(:, ish, :), msh.nqn, msh.nel) .* ...
         u(repmat(sp.connectivity(ish, :), msh.nqn, 1));
     end
     errcurl = errcurl + sum((valnum(:) - curl_valex).^2 .* w);

   case{3}
    for idir = 1:msh.rdim
      valnum = zeros (msh.nqn, msh.nel);
      for ish = 1:sp.nsh_max
        valnum = valnum + ...
          reshape (sp.shape_function_curls(idir, :, ish, :), msh.nqn, msh.nel) .* ...
	      reshape (u(repmat(sp.connectivity(ish,:), msh.nqn, 1)), msh.nqn, msh.nel);
      end
      valex = curl_valex(idir, :);
      errcurl = errcurl + sum ((valnum(:) - valex(:)).^2 .* w);
    end
  end

  errl2  = sp_l2_error (sp, msh, u, uex);
  errhcurl = sqrt (errl2^2 + errcurl);
  errcurl = sqrt (errcurl);

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
clear all; close all;
I=zeros(200, 200);
I(50:150, 50:150)=1;
theta=0:10:180;
[R, xp]=radon(I, theta);
figure;
subplot(121);
imshow(I);
subplot(122);
imagesc(theta, xp, R);
colormap(hot);
colorbar;

author: HuangCongQing


text:
%kvnconddil 'N times (with respect to the composition) the conditional dilation'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vnconddil.pane file
%
% Parameters:
% InputFile: i1 'Input Image #1', required: 'first input image'
% InputFile: i2 'Input Image #2', required: 'second input image'
% InputFile: str 'Struct Element', required: 'Structuring Element'
% Integer: n 'Number of Dilations', default: 5: 'number of cond. dilations'
% OutputFile: o 'Output Image ', required: 'resulting output image'
%
% Example: o = kvnconddil({i1, i2, str}, {'i1','';'i2','';'str','';'n',5;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vnconddil - N times (with respect to the composition) the conditional dilation
%
%  DESCRIPTION
% .I vnconddil
% will do N sucessive cond. dilations of the first input image by a set structuring element, using the second input image as a conditional image.
%
% The two input images must be of the same size and have the same
% number of data bands.  The input
% images may be of the same or different data types.  If the input
% images are of different data types, then all input images will be
% upcast to the highest input image data type.  The output
% image data type will be the same as that of the highest data type
% of the input images.
%
%
%
%  EXAMPLES
% vnconddil -i1 ball.xv -i2 feath.xv -str b272.str -n 5 -o outimage.xv
%
% Will do 5 sucessive cond. dilations of image 1 "ball.xv" by the structuring set b272.str using "feath.xv" as a conditional image, with the resulting image written to "outimage.xv".
%
%  "SEE ALSO"
%
%  RESTRICTIONS
% .I vnconddil
% can be defined for all data types supported by Khoros, but at the moment it has been implemented just for the bit and unsigned char types.
% The structuring elements are subsets of the 3x3 matrix and the origin is always at the center of this matrix.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993-1997 Junior Barrera, Roberto Lotufo.  All rights reserved.
%

function varargout = kvnconddil(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kvnconddil(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'i2', '__input';'str', '__input';'n', 5;'o', '__output'};
maxval={0,0,0,2,0};
minval={0,0,0,2,0};
istoggle=[0,0,0,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','InputFile','Integer','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'vnconddil"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function polygon_properties_test03 ( )

%*****************************************************************************80
%
%% POLYGON_PROPERTIES_TEST03 tests POLYGON_AREA_2;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2014
%
%  Author:
%
%    John Burkardt
%
  test_num = 2;

  area_exact_test = [ 2.0, 6.0 ];
  n_test = [ 4, 8 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYGON_PROPERTIES_TEST03\n' );
  fprintf ( 1, '  For a polygon:\n' );
  fprintf ( 1, '  POLYGON_AREA_2 computes the area.\n' );

  for test = 1 : test_num

    n = n_test(test);
    area_exact = area_exact_test(test);

    if ( test == 1 )

      v = [ ...
        1.0, 0.0; ...
        2.0, 1.0; ...
        1.0, 2.0; ...
        0.0, 1.0 ]';

    elseif ( test == 2 )

      v = [ ...
        0.0, 0.0; ...
        3.0, 0.0; ...
        3.0, 3.0; ...
        2.0, 3.0; ...
        2.0, 1.0; ...
        1.0, 1.0; ...
        1.0, 2.0; ...
        0.0, 2.0 ]';

    end

    r8mat_transpose_print ( 2, n, v, '  The polygon vertices:' );

    area = polygon_area_2 ( n, v );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Exact area is        %g\n', area_exact );
    fprintf ( 1, '  The computed area is %g\n', area );

  end

  return
end


text:
function inv_mix = spm_cfg_eeg_inv_mix
% Configuration file for merging (using a new inversion) a number of
% imaging source inversion reconstructions
%__________________________________________________________________________
% Copyright (C) 2010 Wellcome Trust Centre for Neuroimaging

% Gareth Barnes
% $Id: spm_cfg_eeg_inv_mix.m 5924 2014-03-19 14:59:12Z gareth $

D = cfg_files;
D.tag = 'D';
D.name = 'M/EEG datasets';
D.filter = 'mat';
D.num = [1 Inf];
D.help = {'Select the M/EEG mat files or .mat files containing inverses'};

val = cfg_entry;
val.tag = 'val';
val.name = 'Inversion index';
val.strtype = 'n';
val.help = {'Index of the cell in D.inv (same for all files) where the forward model can be found and the results will be stored.'};
val.val = {1};

prefix = cfg_entry;
prefix.tag = 'prefix';
prefix.name = 'Merged file prefix';
prefix.strtype = 's';
prefix.val={'merged'};
prefix.help = {'Prefix for the new filename that will contain the merged inversion results'};

inv_mix = cfg_exbranch;
inv_mix.tag = 'inv_mix';
inv_mix.name = 'Merge source estimates from multiple inversions';
inv_mix.val = {D, val,prefix};
inv_mix.help = {'To merge different source level variance estimates based on the same data'};
inv_mix.prog = @run_inv_mix;
inv_mix.vout = @vout_inv_mix;
inv_mix.modality = {'MEG'};

function  out = run_inv_mix(job)

inverse = [];
if numel(job.D)<1,
    error('Need to add a number of files to combine');
end;

%% first compile multiple inversions
disp('Loading inversions');
allID=[];
allJ=[];
allqC=[];
allF=[];

surfdir=[];
for j=1:numel(job.D), %% move through files- assume that changing directory means surface(i.e. lead field also changes)
    [a1,b1,c1]=fileparts(job.D{j});
    surfdir=strvcat(surfdir,a1);
end;

gainfiles=[];
for j=1:numel(job.D), %% move through files- assume that changing directory means surface(i.e. lead field also changes)

    try
        spmfilename=job.D{j};
        D = spm_eeg_load(spmfilename);
        inv=D.inv{job.val};
    catch
        dum=load(job.D{j});
        spmfilename=dum.spmfilename;
        [a0,b1,c1]=fileparts(spmfilename);
        [a1,b0,c0]=fileparts(job.D{j});
        D=spm_eeg_load([a1 filesep b1 c1]);
        inv=dum.inv;
    end;

    allF(j)=inv.inverse.F;
    allqC(j,:)=inv.inverse.qC;

    allmesh{j}.M   = inv.mesh.tess_mni;
    allcortex{j}=inv.mesh.tess_ctx;
    %% transform back to temporal subspace of original data
    %% (as Us will be different for different surfaces- but just rotations)
    check_data(j,:)=inv.inverse.U{1}'*inv.inverse.Y(:,1);

    disp('Loading SPM gain matrices from surface directories');

    [L,D]=spm_eeg_lgainmat(D);
    gainfile=[deblank(surfdir(j,:)) filesep D.inv{D.val}.gainmat];
    gainfiles = strvcat(gainfiles,gainfile);
    allL1(j,:)=L(1,:);

end;

%% check original data was the same
if numel(job.D)>1,
    if max(std(check_data))>max(std(check_data'))/1e6,
        error('data is not the same for these files');
    end;
end;

%%NB OCCAMS RAZOR AT 3
useind=find(allF>max(allF)-3); %%
allqC=allqC(useind,:);
allmesh=allmesh(useind);
allcortex=allcortex(useind);
gainfiles=gainfiles(useind,:);
allL1=allL1(useind,:);
%% only take unique lead fields out
[dum,fileind,surfind]=unique(allL1,'rows');
ugainfiles=gainfiles(fileind,:);

[a1,b1,c1]=fileparts(spmfilename);
outdir=surfdir(1,1:max(strfind(surfdir(1,:),filesep))); %% go up one directory
outfilename=[outdir filesep job.prefix b1 c1];
disp(sprintf('Copying and renaming original SPM file to %s',outfilename));

D2=copy(D,outfilename);

if ~isfield(D2.inv{D2.val},'inverse'),
    disp('No inversion parameters in file, taking from last inversion');
    D2.inv{D2.val}.inverse=inv.inverse;
end;

clear dum inverse D

vert=[];
face=[];
cmap1=[];
cortexstr='';
for j=1:length(fileind),
    offset=uint32(size(vert,1).*ones(size(allmesh{fileind(j)}.M.face)));
    col=j*10;
    vert=[vert ;allmesh{fileind(j)}.M.vert];
    face=[face ;allmesh{fileind(j)}.M.face+offset];
    cmap1=[cmap1 ;repmat(col,size(allmesh{fileind(j)}.M.face,1),1)];

    cortexstr=[cortexstr sprintf('%s',allcortex{fileind(j)})];
    if j<length(fileind),
        cortexstr=[cortexstr ';'];
    end;
end;
figure;
h=trisurf(face,vert(:,1),vert(:,2),vert(:,3),cmap1)
set(h,'Linestyle','none');
alpha(0.1);

% UPDATE MNI MESH- MAY NEED TO UPDATE OTHERS TOO
 D2.inv{D2.val}.mesh.tess_mni.vert=vert;
 D2.inv{D2.val}.mesh.tess_mni.face=face;
 D2.inv{D2.val}.mesh.tess_ctx=cortexstr;

%D2.save;
clear vert face cmap1 offset

D2=spm_eeg_invert_classic_mix(D2,D2.val,allqC,surfind,ugainfiles);

disp(sprintf('Improvement in model evidence over best single solution %3.2f',D2.inv{D2.val}.inverse.F-max(allF)));

if ~iscell(D2)
    D2 = {D2};
end

for i = 1:numel(D2)
    save(D2{i});
end

out.D = {outfilename};

function dep = vout_inv_mix(job)
% Output is always in field "D", no matter how job is structured
dep = cfg_dep;
dep.sname = 'M/EEG dataset(s) after imaging source reconstruction';
% reference field "D" from output
dep.src_output = substruct('.','D');
% this can be entered into any evaluated input
dep.tgt_spec   = cfg_findspec({{'filter','mat'}});


text:
function Z=fftfilter(X, H)
F=fft2(X, size(H,1), size(H, 2));
Z=H.*F;
Z=ifftshift(Z);
Z=abs(ifft2(Z));
Z=Z(1:size(X, 1), 1:size(X, 2));
end


text:
function [nInput,inputMin,inputMax, nTarget, targetMin, targetMax] = preNorm(input, target, inputMax)
% Preprocesses data so that minimum is -1 and maximum is 1.
%
%% Syntax
%   [nInput,inputMin,inputMax] = preNorm(input);
%   [nInput,inputMin,inputMax,nTarget,targetMin,targetMax] = preNorm(input, target);
%   [nInput] = preNorm(input,inputMin,inputMax);
%
%
%% Description
% Function normalizes inputs so that they fall in the interval [-1,1].
% Algorithm:  nInput = 2*(input-inputMin)/(inputMax-inputMin) - 1; If the
% data needs to be normalzed using previously computed minimums and
% maximums, we use: nInput = preNorm(input,inputMin,inputMax);
%
% Input:
% * input           ... the n x D input matrix
% * target/inputMin ... the target vector/if there are 3 input arguments it
%                       will be considered as a vector of minimums
% * inputMax        ... maximums
%
% Output:
% * nInput      ... the n x D normalized input matrix
% * inputMin    ... the row vector containing minimums for each dimension
% * inputMax    ... the row vector containing maximums for each dimension
% * nTarget     ... the normalized target vector
% * targetMin   ... the target minimum
% * targetMax   ... the target maximum
%
% See also:
% postNorm, postNormVar

%% Signature
% * Written by Tomaz Sustar, January 2012

if nargin > 3
  error('Wrong number of arguments.');
end

if nargin==3 % normalize targets if given
  inputMin = target;
  [nInput,inputMin,inputMax] = normalize(input,inputMin,inputMax); %normalize input
else
  [nInput,inputMin,inputMax] = normalize(input); %normalize input
end

if nargin==2 % normalize targets if given
  [nTarget, targetMin, targetMax] = normalize(target);
end

function [nInput,inputMin,inputMax] = normalize(input, inputMin, inputMax)
[n, D] = size(input); % n - number of mesurements, D - input space dimenson

if nargin ~= 3
  inputMin = min(input); % row vector of minimiums
  inputMax = max(input); % row vector of maximums
end

isequal = inputMin==inputMax;
notequal = ~isequal;
if sum(isequal) ~= 0
  warning('Some maximums and minimums are equal. Those inputs will not be transformed.');
  inputMin0 = inputMin.*notequal - 1*isequal; % where equal set minimums to -1
  inputMax0 = inputMax.*notequal + 1*isequal; % and maximums to +1 so the data will not be transformed
else
  inputMin0 = inputMin;
  inputMax0 = inputMax;
end

nInput = 2*(input-repmat(inputMin0,n,1))./repmat((inputMax0-inputMin0),n,1) - 1; % normalize

author: Dynamic-Systems-and-GP

repo: GPdyn


text:
function U = fracdiffdemoydelay(alpha,alphad,beta,steps)

% A demo to the article:
% I. Podlubny, A.Chechkin, T. Skovranek, YQ Chen,
% B. M. Vinagre Jara, "Matrix approach to discrete
% fractional calculus II: partial fractional differential
% equations". http://arxiv.org/abs/0811.1355
%
% This function was used in the above article for the
% numerical solution and visualization of Example 5.
% It illustrates the ease of use of the matrix approach
% to discretization of partial differential equations
% with fractional derivatives with respect to the time
% variable with delays. Delays are considered to be multiples
% of a selected time step.
%
% This function solves the time-fractional diffusion-wave equation
% containing one delayed fractional derivative:
%  0.5*u_{t}^(\alpha) (x,t) + 0.5*u_{t}^(\alphad) (x,t-tau*steps)
%        = a2 * u_{xx}^{(2)}(x,t) + f(x,t)
% under zero initial and boundary conditions
%
% For alpha=1, alphad=1, beta=2, steps=0, and f(x,t)=8
% we have the part of the test example from the book:
% W. E. Milne. "Numerical Solution of Differential Equations".
%               New York: Wiley (London: Chapman & Hall), 1953.
%
% which is the classical heat conduction equation
%  u_{t} (x,t) = a2 * u_{xx}(x,t) + f(x,t)

a2=1;       % coefficient from the diffusion equation
L = 1;      % length of spatial interval

% Number of spatial steps + 1 is:
m = 21; % 11, 21

% Number of steps in time + 1 is:
n =148; % 37, 148

h = L / (m-1);          % spatial step
tau = h^2 / (6*a2);     % time step

% generating the matrix for approximation
B1 = ban(alpha,n-1,tau)';       % alpha-th order derivative with respect to time
TD = kron(B1, eye(m));          % time derivative matrix

Bdelay = shift (ban(alphad,n-1+steps,tau)', steps); % delayed derivative of order alphad
TDdelay = kron(Bdelay, eye(m)); % delayed time derivative matrix

B2 = ransym(beta,m,h);          % beta-th order derivative with respect to X
SD = kron(eye(n-1), B2);        % spatial derivative matrix

SystemMatrix = 0.5*TD + 0.5*TDdelay- a2*SD;   % matrix corresponding to discretization
                                      % in space and time

% remove columns with '1' and 'm' from SystemMatrix
S = eliminator (m, [1 m]);
SK = kron(eye(n-1), S);
SystemMatrix_without_columns_1_m = SystemMatrix * SK';

% remove rows with '1' and 'm' from SystemMatrix_without_columns_1_m
S = eliminator (m, [1 m]);
SK = kron(eye(n-1), S);
SystemMatrix_without_rows_columns_1_m = SK * SystemMatrix_without_columns_1_m;

% Right hand side
F = 8*ones(size(SystemMatrix_without_rows_columns_1_m,1),1);

% Solution of the system
Y = SystemMatrix_without_rows_columns_1_m\F;

% Reshape solution array -- values for k-th time step
% are in the k-th column of YS:
YS = reshape(Y,m-2,n-1);
YS = fliplr(YS);

U = YS;

% plot graph
[rows, columns] = size(U);
U = [ zeros(1, columns); U;  zeros(1, columns)];
U = [zeros(1,m)' U];
[XX,YY]=meshgrid(tau*(0:n-1),h*(0:m-1));

mesh(XX,YY,U)
xlabel('t');
ylabel('x');
zlabel('y(x,t)');
title(['\alpha = ', num2str(alpha), ', ', ...
       ' \gamma = ', num2str(alphad), ', ', ...
       ' \beta = ', num2str(beta), ', ', ...
       ' k = ', num2str(steps) ])

set(gca, 'xlim', [0 tau*n], 'zlim', [0 1])
box on


text:
function p = randperm(n)
%RANDPERM Random permutation.
%	RANDPERM(n) is a random permutation of 1:n.
[ignore,p] = sort(rand(1,n));


text:
function [h, compUpVV, compUpVP, compUp] =  lfmComputeH3VV(gamma1_p, gamma1_m, sigma2, t1, ...
    t2, preFactor, mode)

% LFMCOMPUTEH3VV Helper function for computing part of the LFMVXLFMV kernel.
% FORMAT
% DESC computes a portion of the LFMVXLFMV kernel.
% ARG gamma1 : Gamma value for first system.
% ARG gamma2 : Gamma value for second system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG preFactor : precomputed constants.
% ARG mode: indicates the correct precomputations.
% RETURN h : result of this subcomponent of the kernel for the given values.
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

% Evaluation of h

if nargout>1
    [compUpVV{1}, compUpVP{1}, compUp{1}] = lfmvvComputeUpsilonMatrix(gamma1_p,sigma2, t1,t2, mode);
    [compUpVV{2}, compUpVP{2}, compUp{2}] = lfmvvComputeUpsilonMatrix(gamma1_m,sigma2, t1,t2, mode);
    h = preFactor(1)*compUpVV{1} + preFactor(2)*compUpVV{2};
else
    h = preFactor(1)*lfmvvComputeUpsilonMatrix(gamma1_p,sigma2, t1,t2, mode) ...
        + preFactor(2)*lfmvvComputeUpsilonMatrix(gamma1_m,sigma2, t1,t2, mode);

%     %h = preFactor(2)*lfmvvComputeUpsilonMatrix(gamma1_m,sigma2, t1,t2, mode);
%
%      h = lfmvvComputeUpsilonMatrix(gamma1_m,sigma2, t1,t2, mode);
%
%      epsilon = 1e-6;
%      valg = gamma1_m;
%      gamma1_m = valg + epsilon;
%      h1 = lfmvvComputeUpsilonMatrix(gamma1_m,sigma2, t1,t2, mode);
%      gamma1_m = valg - epsilon;
%      h2 = lfmvvComputeUpsilonMatrix(gamma1_m,sigma2, t1,t2, mode);
%      gamma1_m = valg;
%      numerics = 0.5*(h1-h2)/epsilon;
%      theo = lfmvvGradientUpsilonMatrix(gamma1_m,sigma2, t1,t2, mode);

end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% SOLVE_MAXWELL_EIG_MIXED1: Solve the Maxwell eigenvalue problem with a mixed formulation, and a B-spline discretization.
%
% Example to solve the problem
%
%    curl (1/mu(x) curl (u)) = lambda (epsilon(x) u)   in Omega = F((0,1)^n)
%          div (epsilon(x) u) = 0                       in Omega
%       (1/mu(x) curl(u)) x n = 0                       on Gamma_N
%                       u x n = 0                       on Gamma_D
%
% with the variational mixed formulation
%
%    \int (1/mu(x) curl(u) curl(v)) + \int (epsilon(x) v grad(p))
%                = lambda \int (epsilon(x) u v),   \forall v \in H_0(curl),
%                                     \int (epsilon(x) u grad(q)) = 0,
%                                                  \forall q \in H^1_0.
%
% USAGE:
%
%  [geometry, msh, space, sp_mul, eigv, eigf] =
%                  solve_maxwell_eig_mixed1 (problem_data, method_data)
%
% INPUT:
%
%  problem_data: a structure with data of the problem. It contains the fields:
%    - geo_name:     name of the file containing the geometry
%    - nmnn_sides:   sides with Neumann boundary condition (may be empty)
%    - drchlt_sides: sides with Dirichlet boundary condition
%    - c_elec_perm:  electric permittivity (epsilon in the equation)
%    - c_magn_perm:  magnetic permeability (mu in the equation)
%
%  method_data : a structure with discretization data. Its fields are:
%    - degree:     degree of the spline functions.
%    - regularity: continuity of the spline functions.
%    - nsub:       number of subelements with respect to the geometry mesh
%                   (nsub=1 leaves the mesh unchanged)
%    - nquad:      number of points for Gaussian quadrature rule
%
% OUTPUT:
%
%  geometry: geometry structure (see geo_load)
%  msh:      mesh object that defines the quadrature rule (see msh_cartesian)
%  space:    space object that defines the discrete functions (see sp_vector)
%  sp_mul:   space object for the multiplier (see sp_scalar)
%  eigv:     the computed eigenvalues
%  eigf:     degrees of freedom of the associated eigenfunctions
%
% See also EX_MAXWELL_EIG_MIXED1_SQUARE for an example
%
% Copyright (C) 2010, 2011, 2015 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function [geometry, msh, space, sp_mul, eigv, eigf] = ...
              solve_maxwell_eig_mixed1 (problem_data, method_data)

% Extract the fields from the data structures into local variables
data_names = fieldnames (problem_data);
for iopt  = 1:numel (data_names)
  eval ([data_names{iopt} '= problem_data.(data_names{iopt});']);
end
data_names = fieldnames (method_data);
for iopt  = 1:numel (data_names)
  eval ([data_names{iopt} '= method_data.(data_names{iopt});']);
end

% Construct geometry structure
geometry = geo_load (geo_name);

[knots, zeta] = kntrefine (geometry.nurbs.knots, nsub-1, degree, regularity);
[knots_hcurl, degree_hcurl] = knt_derham (knots, degree, 'Hcurl');

% Construct msh structure
rule     = msh_gauss_nodes (nquad);
[qn, qw] = msh_set_quad_nodes (zeta, rule);
msh      = msh_cartesian (zeta, qn, qw, geometry);

% Construct the space structures for the field and the Lagrange multiplier
scalar_spaces = cell (msh.ndim, 1);
for idim = 1:msh.ndim
  scalar_spaces{idim} = sp_bspline (knots_hcurl{idim}, degree_hcurl{idim}, msh);
end
space = sp_vector (scalar_spaces, msh, 'curl-preserving');
sp_mul = sp_bspline (knots, degree, msh);

% Assemble the matrices
if (msh.rdim == 2)
  invmu = @(x,y) 1./c_magn_perm (x,y);
elseif (msh.rdim == 3)
  invmu = @(x,y,z) 1./c_magn_perm (x,y,z);
end

stiff_mat = op_curlu_curlv_tp (space, space, msh, invmu);
mass_mat  = op_u_v_tp (space, space, msh, c_elec_perm);
saddle_mat = op_v_gradp_tp (space, sp_mul, msh, c_elec_perm);

% Apply homogeneous Dirichlet boundary conditions
drchlt_dofs = []; drchlt_dofs_mul = [];
for iside = 1:numel (drchlt_sides)
  drchlt_dofs = union (drchlt_dofs, space.boundary(drchlt_sides(iside)).dofs);
  drchlt_dofs_mul = union (drchlt_dofs_mul, sp_mul.boundary(drchlt_sides(iside)).dofs);
end
if (isempty (drchlt_dofs_mul))
  drchlt_dofs_mul = sp_mul.ndof;
end

int_dofs = setdiff (1:space.ndof, drchlt_dofs);
int_dofs_mul = setdiff (1:sp_mul.ndof, drchlt_dofs_mul);

% Solve the eigenvalue problem
stiff_mat  = stiff_mat (int_dofs, int_dofs);
mass_mat   = mass_mat (int_dofs, int_dofs);
saddle_mat = saddle_mat (int_dofs_mul, int_dofs);

A = [stiff_mat, saddle_mat.'; ...
     saddle_mat, sparse(numel(int_dofs_mul),numel(int_dofs_mul))];
M = [mass_mat, sparse(numel(int_dofs), numel(int_dofs_mul)); ...
     sparse(numel(int_dofs_mul), numel(int_dofs)+numel(int_dofs_mul))];

eigf = zeros (space.ndof + sp_mul.ndof, numel(int_dofs) + numel(int_dofs_mul));
[eigf([int_dofs space.ndof+int_dofs_mul], :), eigv] = eig (full(A), full(M));
eigv = diag (eigv);

end

%!demo
%! ex_maxwell_eig_mixed1_square

%!demo
%! ex_maxwell_eig_mixed1_Lshaped

%!demo
%! ex_maxwell_eig_mixed1_cube

%!demo
%! ex_maxwell_eig_mixed1_thick_L


text:
function sbfmf_write_nframes(fp,version,nframes)

fseek(fp,4+length(version)+4+4,'bof');
fwrite(fp,nframes,'uint32');

author: kristinbranson

repo: JAABA


text:
function bSimilar = isSimilarImage(img1, img2)
%
%
%       bSimilar = isSimilarImage(img1, img2)
%
%
%       This image checks if two images are the same.
%
%     Copyright (C) 2018  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

[r1, c1, ~] = size(img1);
[r2, c2, ~] = size(img2);

bSimilar = (r1 == r2) & (c1 == c2);

end

author: banterle