text:
function varargout = process_import_data_raw( varargin )
% PROCESS_IMPORT_DATA_RAW: Import a raw file in the database.

% @=============================================================================
% 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.
% =============================================================================@
%
% Authors: Francois Tadel, 2012

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Create link to raw file';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = {'Import', 'Import recordings'};
    sProcess.Index       = 11;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/ChannelFile#Link_the_raw_files_to_the_database';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'import'};
    sProcess.OutputTypes = {'raw'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 0;
    sProcess.isSeparator = 1;
    % Option: Subject name
    sProcess.options.subjectname.Comment = 'Subject name:';
    sProcess.options.subjectname.Type    = 'subjectname';
    sProcess.options.subjectname.Value   = 'NewSubject';
    % Option: File to import
    sProcess.options.datafile.Comment = 'File to import:';
    sProcess.options.datafile.Type    = 'datafile';
    sProcess.options.datafile.Value   = {...
        '', ...                                % Filename
        '', ...                                % FileFormat
        'open', ...                            % Dialog type: {open,save}
        'Open raw EEG/MEG recordings...', ...  % Window title
        'ImportData', ...                      % LastUsedDir: {ImportData,ImportChannel,ImportAnat,ExportChannel,ExportData,ExportAnat,ExportProtocol,ExportImage,ExportScript}
        'single', ...                          % Selection mode: {single,multiple}
        'files_and_dirs', ...                  % Selection mode: {files,dirs,files_and_dirs}
        bst_get('FileFilters', 'raw'), ...    % Get all the available file formats
        'DataIn'};                             % DefaultFormats
    % Separator
    sProcess.options.separator.Type = 'separator';
    sProcess.options.separator.Comment = ' ';
    % Replace channel files
    sProcess.options.channelreplace.Comment = 'Replace existing channel file';
    sProcess.options.channelreplace.Type    = 'checkbox';
    sProcess.options.channelreplace.Value   = 1;
    % Align sensors
    sProcess.options.channelalign.Comment = 'Align sensors using headpoints';
    sProcess.options.channelalign.Type    = 'checkbox';
    sProcess.options.channelalign.Value   = 1;
    % Align sensors
    sProcess.options.evtmode.Comment = 'Reading mode (value, bit, ttl, rttl):';
    sProcess.options.evtmode.Type    = 'text';
    sProcess.options.evtmode.Value   = 'value';
    sProcess.options.evtmode.Hidden  = 1;
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    OutputFiles = {};

    % ===== GET OPTIONS =====
    % Get subject name
    SubjectName = file_standardize(sProcess.options.subjectname.Value);
    if isempty(SubjectName)
        bst_report('Error', sProcess, [], 'Subject name is empty.');
        return
    end
    % Get filename to import
    FileName   = sProcess.options.datafile.Value{1};
    FileFormat = sProcess.options.datafile.Value{2};
    if isempty(FileName)
        bst_report('Error', sProcess, [], 'No file selected.');
        return
    end
    % Channel replace
    if sProcess.options.channelreplace.Value
        ChannelReplace = 2;
    else
        ChannelReplace = 0;
    end
    % Channels align
    ChannelAlign = 2 * double(sProcess.options.channelalign.Value);
    % Events reading mode
    if isfield(sProcess.options, 'evtmode') && isfield(sProcess.options.evtmode, 'Value') && ~isempty(sProcess.options.evtmode.Value) && ismember(sProcess.options.evtmode.Value, {'value', 'bit', 'ttl', 'rttl'})
        EventsTrackMode = sProcess.options.evtmode.Value;
    else
        EventsTrackMode = 'value';
    end

    % ===== IMPORT FILES =====
    % Get subject
    [sSubject, iSubject] = bst_get('Subject', SubjectName);
    % Create subject is it does not exist yet
    if isempty(sSubject)
        [sSubject, iSubject] = db_add_subject(SubjectName);
    end
    % Import options
    ImportOptions = db_template('ImportOptions');
    ImportOptions.ChannelReplace  = ChannelReplace;
    ImportOptions.ChannelAlign    = ChannelAlign;
    ImportOptions.DisplayMessages = 0;
    ImportOptions.EventsMode      = 'ignore';
    ImportOptions.EventsTrackMode = EventsTrackMode;
    % Import link to raw
    OutputFiles = import_raw(FileName, FileFormat, iSubject, ImportOptions);
end


text:
% Fig. 5.37  Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
%script to generate figure 5.37, Piper Dakota autopilot design
numG = 160*conv ([1 2.5],[1 0.7]);
denG = conv([1 5 40],[1  .03  .06]);
sysG = tf(numG,denG);
sysD=tf([1   3],[1   20]);
sysDG=sysD*sysG;
figure(1)
rlocus(sysG)
axis([-30 2 -12 12])
grid
hold on
pause;
rlocus(sysDG)
hold off
pause;
figure(2)
rlocus(sysG)
hold on
pause;
rlocus(sysDG)
axis([-1.5 .1    -.6 .6])
grid on
hold off
pause;
K = 0.3;
sysH=tf(1,1);
sysT = feedback (K*sysG,sysH);
sysTD=feedback(sysDG,sysH);
figure(3)
step(sysT)
hold on
step(sysTD)
grid on
hold off

author: Sable

repo: mcbench-benchmarks


text:
classdef EventDeltaVExpendedConstraint < AbstractConstraint
    %EventDeltaVConstraint Summary of this class goes here
    %   Detailed explanation goes here

    properties
        normFact = 1;
        event LaunchVehicleEvent
        eventNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;

        lb(1,1) double = 0;
        ub(1,1) double = 0;

        evalType(1,1) ConstraintEvalTypeEnum = ConstraintEvalTypeEnum.FixedBounds;
        stateCompType(1,1) ConstraintStateComparisonTypeEnum = ConstraintStateComparisonTypeEnum.Equals;
        stateCompEvent LaunchVehicleEvent
        stateCompNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;
    end

    methods
        function obj = EventDeltaVExpendedConstraint(event, lb, ub)
            obj.event = event;
            obj.lb = lb;
            obj.ub = ub;

            obj.id = rand();
        end

        function [lb, ub] = getBounds(obj)
            lb = obj.lb;
            ub = obj.ub;
        end

        function [c, ceq, value, lwrBnd, uprBnd, type, eventNum, valueStateComp] = evalConstraint(obj, stateLog, celBodyData)
            type = obj.getConstraintType();

            deltaVExpended = EventDeltaVExpendedConstraint.computeTotalDeltaV(stateLog, obj.event);
            value = deltaVExpended;

            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                valueStateComp = EventDeltaVExpendedConstraint.computeTotalDeltaV(stateLog, obj.stateCompEvent);

            else
                valueStateComp = NaN;
            end

            [c, ceq] = obj.computeCAndCeqValues(value, valueStateComp);

            lwrBnd = obj.lb;
            uprBnd = obj.ub;

            eventNum = obj.event.getEventNum();
        end

        function sF = getScaleFactor(obj)
            sF = obj.normFact;
        end

        function setScaleFactor(obj, sF)
            obj.normFact = sF;
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesEvent(obj, event)
            tf = obj.event == event;
            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                tf = tf || obj.stateCompEvent == event;
            end
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end

        function tf = usesExtremum(obj, extremum)
            tf = false;
        end

        function tf = canUseSparseOutput(obj)
            tf = false;
        end

        function event = getConstraintEvent(obj)
            event = obj.event;
        end

        function type = getConstraintType(obj)
            type = 'Event Delta-V Expended';
        end

        function [unit, lbLim, ubLim, usesLbUb, usesCelBody, usesRefSc] = getConstraintStaticDetails(obj)
            unit = 'km/s';
            lbLim = 0;
            ubLim = Inf;
            usesLbUb = true;
            usesCelBody = false;
            usesRefSc = false;
        end

        function addConstraintTf = openEditConstraintUI(obj, lvdData)
%             addConstraintTf = lvd_EditGenericMAConstraintGUI(obj, lvdData);

            output = AppDesignerGUIOutput({false});
            lvd_EditGenericMAConstraintGUI_App(obj, lvdData, output);
            addConstraintTf = output.output{1};
        end
    end

    methods(Static, Access=private)
        function deltaVExpended = computeTotalDeltaV(stateLog, event)
            subStateLog = stateLog.getAllStateLogEntriesForEvent(event);

            if(length(subStateLog) > 1)
                g0 = getG0();
                deltaVExpended = 0;

                for(i=1:length(subStateLog)-1) %#ok<NO4LP>
                    stateLogEntry1 = subStateLog(i);
                    stateLogEntry2 = subStateLog(i+1);

                    ut = stateLogEntry1.time;
                    rVect = stateLogEntry1.position;
                    vVect = stateLogEntry1.velocity;

                    bodyInfo = stateLogEntry1.centralBody;
                    tankStates = stateLogEntry1.getAllActiveTankStates();
                    stageStates = stateLogEntry1.stageStates;
                    lvState = stateLogEntry1.lvState;

                    dryMass = stateLogEntry1.getTotalVehicleDryMass();
                    tankStatesMasses = [tankStates.tankMass]';

                    throttleModel = stateLogEntry1.throttleModel;
                    steeringModel = stateLogEntry1.steeringModel;

                    altitude = norm(rVect) - bodyInfo.radius;
                    pressure = getPressureAtAltitude(bodyInfo, altitude);

                    powerStorageStates = stateLogEntry1.getAllActivePwrStorageStates();
                    storageSoCs = NaN(size(powerStorageStates));
                    for(j=1:length(powerStorageStates))
                        storageSoCs(j) = powerStorageStates(j).getStateOfCharge();
                    end

                    throttle = throttleModel.getThrottleAtTime(ut, rVect, vVect, tankStatesMasses, dryMass, stageStates, lvState, tankStates, bodyInfo, storageSoCs, powerStorageStates);

                    [tankMDots, totalThrust, ~] = LaunchVehicleStateLogEntry.getTankMassFlowRatesDueToEngines(tankStates, tankStatesMasses, stageStates, throttle, lvState, pressure, ut, rVect, vVect, bodyInfo, steeringModel, storageSoCs, powerStorageStates, []);

                    if(abs(sum(tankMDots)) > 0)
                        tankMDotsKgS = tankMDots * 1000;
                        totalMDotKgS = sum(tankMDotsKgS); %should be negative
                        totalThrustN = totalThrust * 1000;
                        effIsp = totalThrustN / (getG0() * abs(totalMDotKgS)); %sec

                        totalMass1 = dryMass + stateLogEntry1.getTotalVehiclePropMass();
                        totalMass2 = dryMass + stateLogEntry2.getTotalVehiclePropMass();

                        if(totalMass1 > totalMass2 && stateLogEntry1.time < stateLogEntry2.time)
                            deltaVExpended = deltaVExpended + (g0 * effIsp * log(totalMass1 / totalMass2))/1000; %need to convert to km/s
                        end
                    end
                end
            else
                deltaVExpended = 0;
            end
        end
    end

    methods(Static)
        function constraint = getDefaultConstraint(~, ~)
            constraint = EventDeltaVExpendedConstraint(LaunchVehicleEvent.empty(1,0),0,0);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% Use get_pencil_curves to draw a polyline for an elaborate two-dimensional
% shape that you like. Be creative!

author: odedstein

repo: sgi-introduction-course


text:
function [U S] = performPCA(X)

% This function calculates the covariance matrix (Sigma) and then performs
% svd on it to get the eigen vectors and eigen values

[m n] = size(X); % m - no of egs; n - no of features

Sigma = (1/m)*(X')*(X); % Covariance matrix

[U S V] = svd(Sigma); % Perform Singular Value Decomposition


text:
function im = drawlineonimage(im,x,y,color)

if ~exist('color'),
  color = lines(1);
end

if max(im(:)) > 1,
  color = color*255;
end

npts = length(x);
if ndims(im) < 3,
  im = repmat(im,[1,1,3]);
end

[nr,nc,three] = size(im);
im = mat2cell(im,nr,nc,ones(1,3));

for i = 1:npts-1,

  dx = x(i+1) - x(i);
  dy = y(i+1) - y(i);
  m = dy / dx;
  dx = abs(dx); dy = abs(dy);
  if dx > dy,
    xcurr = linspace(x(i),x(i+1),dx);
    ycurr = y(i) + m*(xcurr-x(i));
  else,
    ycurr = linspace(y(i),y(i+1),dy);
    xcurr = x(i) + (ycurr - y(i))/m;
  end
  for rx = 1:2,
    if rx == 1,
      roundx = floor(xcurr);
    else,
      roundx = ceil(xcurr);
    end
    dx = xcurr - roundx;
    for ry = 1:2,
      if ry == 1,
        roundy = floor(ycurr);
      else,
        roundy = ceil(ycurr);
      end
      dy = ycurr - roundy;
      w = max(0,sqrt(dx.^2 + dy.^2) - (sqrt(2) - 1));
      for c = 1:3,
        im{c}(sub2ind([nr,nc],roundy,roundx)) = im{c}(sub2ind([nr,nc],roundy,roundx)).*w + ...
            color(c)*(1-w);
      end
    end
  end
end

im = cell2mat(im);

author: kristinbranson

repo: JAABA


text:
function hermite_polynomial_test07 ( )

%*****************************************************************************80
%
%% HERMITE_POLYNOMIAL_TEST07 tests HE_QUADRATURE_RULE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 March 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_POLYNOMIAL_TEST07:\n' );
  fprintf ( 1, '  HE_QUADRATURE_RULE computes the quadrature rule\n' );
  fprintf ( 1, '  associated with He(n,x);\n' );

  n = 7;
  [ x, w ] = he_quadrature_rule ( n );

  r8vec2_print ( n, x, w, '      X            W' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Use the quadrature rule to estimate:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    Q = Integral ( -oo < X < +00 ) X^E exp(-X^2) dx\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   E       Q_Estimate      Q_Exact\n' );
  fprintf ( 1, '\n' );

  for e = 0 : 2 * n - 1
    if ( e == 0 )
      f = ones ( n, 1 );
    else
      f(1:n) = x(1:n).^e;
    end
    q = w' * f;
    q_exact = he_integral ( e );
    fprintf ( 1, '  %2d  %14g  %14g\n', e, q, q_exact );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% RPCA | MoG-RPCA | Mixture of Gaussians RPCA (Zhao et al. 2014)
% process_video('RPCA', 'MoG-RPCA', 'dataset/demo.avi', 'output/demo_MoG-RPCA.avi');

%{
clear, clc;
load('dataset/trafficdb/traffic_patches.mat');
V = im2double(imgdb{100});
[M,m,n,p] = convert_video3d_to_2d(V);
%}

r = 1;
param.mog_k = 3;
param.lr_init = 'SVD';
param.maxiter = 100;
param.initial_rank = 2*r;
param.tol = 1e-3;

lr_prior.a0 = 1e-6;
lr_prior.b0 = 1e-6;

mog_prior.mu0 = 0;
mog_prior.c0 = 1e-3;
mog_prior.d0 = 1e-3;
mog_prior.alpha0 = 1e-3;
mog_prior.beta0 = 1e-3;

[lr_model, mog_model, r] = mog_rpca(M, param, lr_prior, mog_prior);

L = lr_model.U*lr_model.V';
S = M - L;

%{
show_2dvideo(M,m,n);
show_2dvideo(L,m,n);
show_2dvideo(S,m,n);
show_2dvideo(hard_threshold(S),m,n);
%}

author: andrewssobral

repo: lrslibrary


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% REAL-TIME BCILAB/SIFT INTERFACE DEMO SCRIPT
%
% This script demonstrates a beta interface between BCILAB
% (Christian Kothe: www.sccn.ucsd.edu/wiki/BCILAB) and SIFT (Tim Mullen:
% www.sccn.ucsd.edu/wiki/SIFT) to perform real-time data extraction,
% cleaning, connectivity modeling and visualization.
%
% This script was prepared for the following whitepaper/demonstration:
% [1] Tim Mullen, Christian Kothe, Yu Chi, Tzyy-Ping Jung,
%     and Scott Makeig. "Real-Time Estimation and 3D Visualization of Sparse
%     Multivariate Effective Connectivity Using Wearable EEG." July 15th,
%     IEEE EMB/CAS/SMC Workshop on Brain-Machine-Body Interfaces. 34th
%     Annual International IEEE EMBS Conference of the IEEE Engineering in
%     Medicine and Biology Society (IEEE EMBS).
%
% Please refer to [1] when using any functionality from this demo.
%
% Copyright Tim Mullen, Christian Kothe, July 2012, SCCN/INC/UCSD
%
% Author: Tim Mullen,       July 2012, SCCN/INC/UCSD
%         Christian Kothe,  July 2012, SCCN/INC/UCSD
%
% This function is part of the Source Information Flow Toolbox (SIFT)
%
% 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% SET UP CONFIGURATION OPTIONS
RUN_LSL          = false;           % If RUN_LSL = true, then stream 'online' from device; If RUN_LSL=false then playback TestingDataFile (below)

% Source reconstruction options (leave disabled)
COLOR_SOURCE_ROI = true;          % this will use special meshes for coloring ROIs
HEAD_MODEL_NAME  = 'resources:/headmodels/Cognionics_64_HeadModelObj_3751.mat';             % path to head model object for source reconstruction. Leave empty if you aren't doing source reconstruction

% Establish file paths
% NOTE: all paths are relative to 'datapath' which is a
% platform-independent path which itself can be relative to bcilab root
% folder (i.e. data:/ is the userdata folder in the bcilab root dir)
datapath         = 'data:/';       % this is relative to the BCILAB root dir
TrainingDataFile = 'Cognionics_64_training.set'; %'calibration_mindo.xdf'; % %'calibration.xdf'; %'Cognionics_Pyramind_demo.set'; %'clean_reversed.xdf'; %'noisy.xdf'; %'Cognionics_Pyramind_demo.set';             % this is the relative path to the calibration dataset
TestingDataFile  = 'Cognionics_64_testing.set'; %'testing.xdf'; %'Cognionics_Pyramind_demo.set'; %'clean_reversed.xdf'; %'noisy.xdf'; %'Cognionics_Pyramind_demo.set';             % this is an optional path to a dataset to playback (if RUN_LSL = false)
GUI_CONFIG_NAME  = 'DEMO_SOURCELOC_METACP_CFG.mat'; %'DARPA_DEMO_METACP_CFG_FEWCHANS.mat';             % relative path to a default pipeline configuration
GUI_BRAINMOVIE_CONFIG_NAME = 'DEMO_SOURCELOC_BM_CFG.mat'; %'DARPA_DEMO_BM_CFG.mat';   % relative path to BrainMovie configuration

% Set up the name of the stream we will write to in the workspace
streamname              = 'EEGDATA';
LSL_SelectionProperty   = 'type';  % can also be 'name'
LSL_SelectionValue      = 'EEG';

%% load head model object
if ~isempty(HEAD_MODEL_NAME)
    hmObj = headModel.loadFromFile(env_translatepath(HEAD_MODEL_NAME));
    % load meshes for visualization
    surfData    = load(hmObj.surfaces);
    surfData    = surfData.surfData;
    scalpMesh   = surfData(1);
    csfMesh     = surfData(2);
    cortexMesh  = surfData(3);
end

%% load calibration recording
calibData = exp_eval(io_loadset([datapath TrainingDataFile],'markerchannel',{'remove_eventchns',false}));
flt_pipeline('update');

%% start LSL streaming
if RUN_LSL == true
    run_readlsl('MatlabStream',streamname,'SelectionProperty',LSL_SelectionProperty,'SelectionValue',LSL_SelectionValue);
end
%% ... OR read from datafile
if RUN_LSL == false
    run_readdataset('MatlabStream',streamname,'Dataset',io_loadset([datapath TestingDataFile],'markerchannel',{'remove_eventchns',false}));
end

%% initialize some variables
[benchmarking.preproc     ...
 benchmarking.modeling    ...
 benchmarking.brainmovie] = deal(NaN);

newPipeline     = true;
newBrainmovie   = true;
specOpts        = [];

% create a minimal moving average pipeline...
% noflats = exp_eval_optimized(flt_movavg(flt_selchans(calibData,{'FC3', 'FC4', 'FC2', 'FC1', 'C3', 'C4', 'C2', 'C1', 'CP3', 'CP4', 'CP2', 'CP1', 'P3', 'P4', 'P2', 'P1'})));
noflats = flt_seltypes(calibData,'EEG');
%noflats = exp_eval_optimized(flt_movavg(flt_seltypes(calibData,'EEG')));
highpassed   = exp_eval_optimized(flt_fir(noflats,[0.5 1],'highpass'));
goodchannels = exp_eval_optimized(flt_clean_channels('signal',highpassed,'min_corr',0.35,'ignored_quantile',0.3));

%% initialize options
if ~isempty(GUI_CONFIG_NAME) && exist(env_translatepath([datapath GUI_CONFIG_NAME]),'file')
    % try to load a configuration file
    io_load([datapath GUI_CONFIG_NAME]);
else
    % manually set some defaults
    opts.miscOptCfg.dispBenchmark  = true;
    opts.miscOptCfg.doSIFT         = [];
    opts.holdPipeline              = false;
    opts.exitPipeline              = false;

    opts.siftPipCfg.preproc.verb = 0;

    % set some defaults
    opts.fltPipCfg = ...
        struct('DataCleaning', ...
                struct('RetainPhases',true, ...
                       'HaveBursts', true, ...
                       'HaveChannelDropouts', false, ...
                       'DataSetting', 'drycap'), ...
                'Rereferencing', {{}});

    opts.siftPipCfg = ...
        struct('preproc', ...
                    struct('normalize', ...
                        struct('verb', 0, ...
                               'method', {{'time'}})));

end

% additional options regarding brainmovie node/edge adaptation
opts.adaptation.adaptationHL    = 10;  % HL of exp. win. MA in frames
opts.adaptation.updateInterval  = 1;
opts.adaptation.bufferTime      = 1;

% set the arg_direct flag to true
opts = arg_setdirect(opts,true);

%% initialize the Meta Control Panel
if ~isempty(HEAD_MODEL_NAME)
    calibData.srcpot = 1; % allow sources
end
figHandles.MetaControlPanel = gui_metaControlPanel(calibData,opts);
waitfor(figHandles.MetaControlPanel,'UserData','initialized');

%% Initialize figures
figHandles.specDisplay      = [];
figHandles.BMDisplay        = [];
figHandles.BMControlPanel   = [];

%% try to load a BrainMovie configuration
try
    io_load([datapath GUI_BRAINMOVIE_CONFIG_NAME]);
    BMCFG.BMopts.bmopts_suppl = {'title',{'Multivariate ','Granger Causality  '}};
    BMCFG.BMopts.caption = true;
catch
    BMCFG = struct([]);
end

BMRenderMode = 'init_and_render';

%% Main loop
% -------------------------------------------------------------------------
while ~opts.exitPipeline

    try

        if newPipeline
            % create a new pipeline on training data
            fprintf('Pipeline changed\n');
            cleaned_data = exp_eval(flt_pipeline('signal',calibData,opts.fltPipCfg));
            pipeline     = onl_newpipeline(cleaned_data,streamname);
            newPipeline  = false;

%             disp('profile on');
%             profile -memory on
        end

        if opts.holdPipeline
            % pause the pipeline
            pause(1);
            continue;
        end

        % grab a chunk of data from the stream and preprocess it
        % ---------------------------------------------------------------------
        prepbench = tic;
        [eeg_chunk,pipeline] = onl_filtered(pipeline, ...
                                            round(opts.miscOptCfg.winLenSec*cleaned_data.srate), ...
                                            opts.miscOptCfg.suppress_console_output, ...
                                            fastif(isfield(opts.fltPipCfg,'psourceLocalize') && opts.fltPipCfg.psourceLocalize.arg_selection,{'srcpot'},{}));
        benchmarking.preproc = toc(prepbench);

        % model the chunk via sift
        % ---------------------------------------------------------------------
        if opts.miscOptCfg.doSIFT.arg_selection
            try
                % force the modeling window length to agree with chunk length
                opts.siftPipCfg.modeling.winlen = opts.miscOptCfg.winLenSec;

                % import ROI Names
                if strcmpi(opts.siftPipCfg.preproc.sigtype.arg_selection,'sources')
                    opts.siftPipCfg.preproc.varnames = eeg_chunk.roiLabels;
                end

                % select a subset of channels
                if ~isempty(opts.miscOptCfg.doSIFT.channelSubset)
                    [dummy eeg_chunk] = evalc('hlp_scope({''disable_expressions'',1},@flt_selchans,''signal'',eeg_chunk,''channels'',opts.miscOptCfg.doSIFT.channelSubset,''arg_direct'',true)');
                end

                % run the SIFT modeling pipeline
                siftbench = tic;
                [eeg_chunk cfg] = hlp_scope({'is_online',1},@onl_siftpipeline,'EEG',eeg_chunk,opts.siftPipCfg);
                benchmarking.modeling = toc(siftbench);

                if eeg_chunk.pnts == 0
                    fprintf('No data!\n');
                    pause(0.01);
                    continue;
                end

            catch err
                hlp_handleerror(err);
                pause(0.1);
                continue;
            end

            % visualize the brainmovie
            % -------------------------------------------------------------
            if opts.miscOptCfg.doSIFT.dispBrainMovie.arg_selection

                if newBrainmovie ...
                    | isempty(ishandle(figHandles.BMDisplay)) ...
                    | ~ishandle(figHandles.BMDisplay)

                    % initialize brainmovie
                    if ~isempty(BMCFG)
                        BMCFG.showNodeLabels.nodelabels = eeg_chunk.CAT.curComponentNames;
                        figHandles.BMControlPanel = gui_causalBrainMovie3D_online(eeg_chunk,eeg_chunk.CAT.Conn,struct('arg_direct',0),BMCFG,'timeRange',[]);
                    else
                        figHandles.BMControlPanel = gui_causalBrainMovie3D_online(eeg_chunk,eeg_chunk.CAT.Conn,struct('arg_direct',0),'showNodeLabels',{'nodelabels',eeg_chunk.CAT.curComponentNames});
                    end
                    waitfor(figHandles.BMControlPanel,'UserData','newBrainmovie');

                    % set the closing behavior
                    set(figHandles.BMControlPanel,'CloseRequestFcn','evalin(''base'',''figHandles.BMControlPanel = [];''); delete(gcbf);');

                    BM_CFG_CHANGED       = true;
                    newBrainmovie        = false;
                    figHandles.BMDisplay = [];
                    BMStateVars          = [];
                end

                if BM_CFG_CHANGED
                    BMRenderMode    = 'init_and_render';
                    BMStateVars     = [];

                    % determine whether to adapt limits
                    adapt_nodeSizeDataRange  = isempty(BMCFG.BMopts.graphColorAndScaling.nodeSizeDataRange)  && opts.miscOptCfg.doSIFT.dispBrainMovie.adaptBrainMovieLimits;
                    adapt_edgeSizeDataRange  = isempty(BMCFG.BMopts.graphColorAndScaling.edgeSizeDataRange)  && opts.miscOptCfg.doSIFT.dispBrainMovie.adaptBrainMovieLimits;
                    adapt_nodeColorDataRange = isempty(BMCFG.BMopts.graphColorAndScaling.nodeColorDataRange) && opts.miscOptCfg.doSIFT.dispBrainMovie.adaptBrainMovieLimits;
                    adapt_edgeColorDataRange = isempty(BMCFG.BMopts.graphColorAndScaling.edgeColorDataRange) && opts.miscOptCfg.doSIFT.dispBrainMovie.adaptBrainMovieLimits;

                    % clear state
                    nodeSizeState   = [];
                    nodeColorState  = [];
                    edgeSizeState   = [];
                    edgeColorState  = [];

                    % reset flag
                    BM_CFG_CHANGED = false;
                end

                % Handle special rendering of source meshes and colors
                % ---------------------------------------------------------
                if COLOR_SOURCE_ROI
                    BG_COLOR = [0.1 0.1 0.1];

                    if BMCFG.BMopts.Layers.custom.arg_selection
                        BMCFG.BMopts.Layers.custom.volumefile = eeg_chunk.dipfit.surfmesh;
                        BMCFG.BMopts.Layers.custom.meshcolor  = hlp_getROIVertexColorTable( ...
                                size(eeg_chunk.dipfit.surfmesh.vertices,1), ...
                                eeg_chunk.roiVertices,BG_COLOR,@(x)distinguishable_colors(x,BG_COLOR));
                    end

                    % if the cortex mesh is a custom mesh with 'constant'
                    % coloring...
                    if strcmpi(BMCFG.BMopts.Layers.cortex.cortexres,'custommesh') ...
                       && strcmp(BMCFG.BMopts.Layers.cortex.cortexcolor.arg_selection,'Constant')
                            % ... replace the colors with ROI coloring
                            BMCFG.BMopts.Layers.cortex.cortexcolor.colormapping = ...
                                hlp_getROIVertexColorTable( ...
                                    size(eeg_chunk.dipfit.surfmesh.vertices,1), ...
                                    eeg_chunk.roiVertices,BG_COLOR,@(x)distinguishable_colors(x,BG_COLOR));
%                         {hmObj.atlas.label, hlp_getROIColorTable(hmObj.atlas.label,eeg_chunk.roiLabels,[0.5 0.5 0.5],[1 0 0])};
                    end
                end

                % Render the Brain Movie
                % ---------------------------------------------------------
                bmbench = tic;
                [BMcfg_tmp BMhandles BMout] = ...
                    vis_causalBrainMovie3D(eeg_chunk,eeg_chunk.CAT.Conn,BMCFG, ...
                        'timeRange',[], ...
                        'BMopts',hlp_mergeVarargin(BMCFG.BMopts, ...
                                    'figurehandle',figHandles.BMDisplay, ...
                                    'RenderMode',BMRenderMode, ...
                                    'speedy',true, ...
                                    'InternalStateVariables',BMStateVars));
                benchmarking.brainmovie = toc(bmbench);

                if isempty(BMStateVars)
                    set(figHandles.BMDisplay, 'MenuBar','figure', ...
                                'CloseRequestFcn','evalin(''base'',''figHandles.BMDisplay = [];''); delete(gcbf);', ...
                                'ToolBar','none', 'Name','BrainMovie3D');
                end

                % save figure handle and state
                figHandles.BMDisplay = BMhandles.figurehandle;
                BMStateVars          = BMcfg_tmp.BMopts.vars;
                BMRenderMode         = 'render';

                % adapt the node and edge color/size limits
                % ---------------------------------------------------------
                if opts.miscOptCfg.doSIFT.dispBrainMovie.arg_selection ...
                   && opts.miscOptCfg.doSIFT.dispBrainMovie.adaptBrainMovieLimits

                    % adapt BrainMovie node/edge limits
                    if adapt_nodeSizeDataRange
                        [BMCFG.BMopts.graphColorAndScaling.nodeSizeDataRange, nodeSizeState] ...
                            = hlp_adaptMinMaxLimits('values',   BMout.NodeSize,     ...
                                                    'state',    nodeSizeState,      ...
                                                    'adaptOpts',opts.adaptation);
                    end
                    if adapt_nodeColorDataRange
                        [BMCFG.BMopts.graphColorAndScaling.nodeColorDataRange, nodeColorState] ...
                            = hlp_adaptMinMaxLimits('values',   BMout.NodeColor,     ...
                                                    'state',    nodeColorState,      ...
                                                    'adaptOpts',opts.adaptation);
                    end
                    if adapt_edgeSizeDataRange
                        [BMCFG.BMopts.graphColorAndScaling.edgeSizeDataRange, edgeSizeState] ...
                            = hlp_adaptMinMaxLimits('values',   BMout.EdgeSize,     ...
                                                    'state',    edgeSizeState,      ...
                                                    'adaptOpts',opts.adaptation);
                    end
                    if adapt_edgeColorDataRange
                        [BMCFG.BMopts.graphColorAndScaling.edgeColorDataRange, edgeColorState] ...
                            = hlp_adaptMinMaxLimits('values',   BMout.EdgeColor,     ...
                                                    'state',    edgeColorState,      ...
                                                    'adaptOpts',opts.adaptation);
                    end

                end % adaptBrainMovieLimits

            else % dispBrainMovie == false
                % close the brainmovie figure if control panel is closed
                if isempty(figHandles.BMControlPanel) ...
                        || ~ishandle(figHandles.BMControlPanel) ...
                        && ishandle(figHandles.BMDisplay)
                    close(figHandles.BMDisplay);
                end
                benchmarking.brainmovie = NaN;
            end

            % display the VAR spectrum
            % -------------------------------------------------------------
            if opts.miscOptCfg.doSIFT.dispSpectrum
                try
                    if ishandle(figHandles.specDisplay)
                        % update the SpecViewer
                        vis_autospectrum(specOpts,'FigureHandle',figHandles.specDisplay,'stream',eeg_chunk);
                    else
                        % generate a GUI and initialize SpecViewer
                        [figHandles.specDisplay specOpts] = arg_guidialog(@vis_autospectrum,'Parameters',{'stream',eeg_chunk});

                        set(figHandles.specDisplay,'CloseRequestFcn','evalin(''base'',''opts.miscOptCfg.dispSpectrum=false; figHandles.specDisplay = [];''); delete(gcbf);');
                        specOpts.arg_direct = 0;
                    end
                catch e
                    disp(e.message);
                    opts.miscOptCfg.doSIFT.dispSpectrum = false;
                end
            else
                if ishandle(figHandles.specDisplay)
                    close(figHandles.specDisplay);
                end
            end

        else
            benchmarking.modeling = NaN;
        end

        if opts.miscOptCfg.dispBenchmark
            % TODO: implement this as a figure
            vis_benchmark(benchmarking);
        end

    catch err
         hlp_handleerror(err);
    end

    pause(0.005);

end

%% all done, close all open figures
fnames = fieldnames(figHandles);
for f=1:length(fnames)
    if ishandle(figHandles.(fnames{f}))
        close(figHandles.(fnames{f}));
    end
end

repo: REST


text:
function [ y, m, d, h, n, s, ierror ] = ymdhms_check_common ( y, m, d, h, n, s )

%*****************************************************************************80
%
%% YMDHMS_CHECK_COMMON checks a Common YMDHMS date.
%
%  Discussion:
%
%    The routine will correct certain simple errors in dates, such as
%      "11:03:42 31 September 1996"
%    which will become
%      "11:03:42 1 October 1996".
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, M, D, H, N, S.
%    These items have the obvious meanings.
%    The routine may change any of these values to more reasonable values.
%
%    Output, integer IERROR, is 0 if no error was detected in the
%    date, and 1 otherwise.
%

%
%  Check that the second is between 0 and 59.
%  N may get bumped up or down.
%
  [ y, m, d, h, n, s ] = second_borrow_common ( y, m, d, h, n, s );

  [ y, m, d, h, n, s ] = second_carry_common ( y, m, d, h, n, s );
%
%  Check that the minute is between 0 and 59.
%  H may get bumped up or down.
%
  [ y, m, d, h, n ] = minute_borrow_common ( y, m, d, h, n );

  [ y, m, d, h, n ] = minute_carry_common ( y, m, d, h, n );
%
%  Check that the hour is between 0 and 23.
%  D may get bumped up or down.
%
  [ y, m, d, h ] = hour_borrow_common ( y, m, d, h );

  [ y, m, d, h ] = hour_carry_common ( y, m, d, h );
%
%  Now make adjustments to D, M, and Y.
%
  [ y, m, d, ierror ] = ymd_check_common ( y, m, d );

  return
end


text:
 function displace = penalty_displace(offsets, sizes)
%function displace = penalty_displace(offsets, sizes)
%|
%| Convert scalar offsets to vector displacements, i.e.,
%| find d1,d2,... such that d1*1 + d2*n1 + d3*n1*n2 + ... = offset.
%| Example: if offset = n1, then [d1 d2 d3] = [0 1 0]
%| Assumes that d <= floor(n/2) for n > 2
%|
%| in
%|	offsets	[LL]		see penalty_offsets.m
%|	sizes	[1 ndim]
%|
%| out
%|	displace [LL ndim]	[d1 d2 ...] for each offset
%|
%| Copyright 2006-12-6, Jeff Fessler, University of Michigan

if nargin == 1 && streq(offsets, 'test'), penalty_displace_test, return, end
if nargin < 2, ir_usage, end

ndim = numel(sizes);
offsets = offsets(:); % [LL 1]
displace = nan(length(offsets), ndim);

if 1 % simple way that assumes d / n < 1/2, which is fine for n > 2
	if length(sizes) == 2 && sizes(1) == 2 % [2 ny] special case!
		if isequal(offsets, [1 2 3 1]') % [1 nx nx+1 nx-1]
			displace = [1 0; 0 1; 1 1; -1 1];
		elseif isequal(offsets, [1]) % [1]
			displace = [1 0];
			warn 'caution: offset=1 is ambiguous when n1=2'
		elseif isequal(offsets, [2]) % [n1]
			displace = [0 1];
		elseif isequal(offsets, [1 2]') % [2 n2]
			displace = [1 0; 0 1];
		else
			pr offsets
			fail 'not done'
		end

	elseif any(sizes(1:end-1) == 2)
		fail 'not done: ambiguous'

	else
		residual = offsets;
		for id=ndim:-1:1
			nd = prod(sizes(1:id-1));
			displace(:,id) = round(residual / nd);
			residual = residual - displace(:,id) * nd;
		end
		if any(residual ~= 0), fail 'bug', end
	end

else

	half = max(floor(sizes/2), 1); % upper bound on dx is almost half the size
	%half = sizes/2;

	for ll=1:length(offsets)
		offset = double(offsets(ll)); % trick: necessary!
		subval = 0;
		dis = zeros(1,ndim);
		for id=ndim:-1:2
			tmp = offset + sum((half(1:id-1)-1) .* [1 sizes(1:id-2)]);
			tmp = tmp + prod(sizes(1:id)) - subval;
			dis(id) = floor(tmp / prod(sizes(1:id-1))) - sizes(id);
			subval = subval + dis(id) * prod(sizes(1:id-1));
		end
		dis(1) = offset - subval;

%		if ~(all(abs(dis) < half || (dis == 1 && half == 1)))
		if ~all(abs(dis) <= half)
			pr dis
			pr half
			warn 'bug?'
		end
		displace(ll,:) = dis;
	end
end

offset_check = displace * [1 cumprod(sizes(1:end-1))]';
jf_equal(offsets, offset_check)

end % penalty_displace

% penalty_displace_test()
function penalty_displace_test

jf_equal(penalty_displace([], [1]), nan(0,1))
jf_equal(penalty_displace([], [1 1]), nan(0,2))

if 1 % check tiny
	% when nx=ny=2 then [1 nx nx+1 nx-1] = [1 2 3 1] so it is ambiguous!
	nx = 2; ny = 2;
	offsets = [1 nx nx+1 nx-1];
	dd = penalty_displace(offsets, [nx ny]);
	jf_equal(dd, [1 0; 0 1; 1 1; -1 1])

	nx = 2; ny = 3;
	dd = penalty_displace(offsets, [nx ny]);
	jf_equal(dd, [1 0; 0 1; 1 1; -1 1])
end

nx = 3; ny = 3;
offsets = [1 nx nx+1 nx-1];
dd = penalty_displace(offsets, [nx ny]);
jf_equal(dd, [1 0; 0 1; 1 1; -1 1])

nx = 100; ny = 80; % 2d
[ix iy] = ndgrid(-2:2, -2:2);
offsets = col(ix + iy * nx);
dd = penalty_displace(offsets, [nx ny]);
jf_equal(dd, [ix(:) iy(:)])

nx = 10; ny = 8; nz = 7; % 3d
[ix iy iz] = ndgrid(-2:2, -2:2, -2:2);
offsets = col(ix + iy * nx + iz * nx * ny);
dd = penalty_displace(offsets, [nx ny nz]);
jf_equal(dd, [ix(:) iy(:) iz(:)])

dim = [8 9]; % for won
offsets2 = penalty_offsets('2d:hvd', dim);
dd2 = penalty_displace(offsets2, dim);

for nz=1:3
	dim3 = [dim nz];
	dd3 = penalty_displace(offsets2, dim3);
	jf_equal(dd3, [dd2 zeros(4,1)])
end

if 1
	dim3 = [10 20 30];
	offsets = penalty_offsets('3d:26', dim3);
	[[1:13]', penalty_displace(offsets, dim3)]
end

end % penalty_displace_test()

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function fig_handle = montage(image_obj, varargin)
% Create a montage of an image_vector (or statistic_image or fmri_data)
% object, on top of a standard anatomical underlay.
% - Takes optional inputs to fmridisplay.addblobs (see help for options)
% - Max 4 images in object. Use get_wh_image to select images from object if needed.
%
% *Usage:
% ::
%
%    [fig_handle or o2 fmridisp object] = montage(image_obj, [optional arguments])
%
% :Optional inputs:
%   **fmridisplay:**
%        for fmridisplay object style montage [default]
%
%   **scnmontage:**
%        for circa 2008-style SCN lab montage for each image vector
%
% :Examples:
% ::
%
%    o2 = montage(mask);
%
%    o2 = montage(dat, 'trans');
%    o2 = montage(dat, 'trans', 'color', [1 0 0]);
%    o2 = montage(dat, 'trans', 'color', [1 0 0], 'transvalue', .4);
%    o2 = montage(dat, 'trans', 'maxcolor', [1 .3 0], 'mincolor', [.5 0 1]);
%    o2 = montage(dat, 'trans', 'mincolor', [.5 0 1], 'transvalue', .5);
%
%    Options for canlab_results_fmridisplay are also passed forward, so that
%    different named montage types can be used.  E.g.,:
%
%    o2 = montage(t, 'trans', 'full');
%
% Set all color maps to the same range:
% o2 = montage(dat, 'trans', 'mincolor', [.5 0 1], 'transvalue', .7, 'cmaprange', [0 3]);

meth = 'fmridisplay';
% montagetype = 'compact2';  % default. This is passed into canlab_results_fmridisplay

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}
            case 'scnmontage', meth = 'scnmontage';

            case {'trans', 'color' 'maxcolor', 'mincolor', 'transvalue', 'cmaprange', 'full', 'compact2' 'noverbose', 'indexmap'}
                % ignore these - passed through

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

% number of images (cols)
n = size(image_obj.dat, 2);

switch meth
    case 'fmridisplay'

        if n > 4
            disp('Warning: Showing first 4 images in data object only.');
            n = 4;
        end

        % r = region(image_obj);

        if n > 1
            do_multirow = true;
            o2 = canlab_results_fmridisplay([], 'noblobs', 'nooutline', 'multirow', n, varargin{:});

            % scale so all maps are on same color scale
            [~, cmaprange] = get_data_range(image_obj);

        else
            do_multirow = false;
            o2 = canlab_results_fmridisplay([], 'noblobs', 'nooutline', varargin{:}); % pass forward args.
        end

        montage_indx = 1:2:2*n;

        for i = 1:n

            obj = get_wh_image(image_obj, i);

            if do_multirow % plot only on montages for this image

                % add title - new in 2021, so use try...catch for now
                if isa(image_obj, 'statistic_image') && ~isempty(image_obj.image_labels)
                    try
                        o2 = title_montage(o2, montage_indx(i)+1, image_obj.image_labels{i});
                    catch
                        disp('Problem displaying title montage for statistic_image object');
                    end
                end

                o2 = addblobs(o2, region(obj, 'noverbose'), 'cmaprange', cmaprange, 'nooutline', varargin{:}, 'wh_montages', [montage_indx(i) montage_indx(i)+1]);

            else % just one image, plot on all montages

                % add title - new in 2021, so use try...catch for now
                if isa(image_obj, 'statistic_image') && ~isempty(image_obj.image_labels)
                    try
                        o2 = title_montage(o2, 5, image_obj.image_labels{1});
                    catch
                        disp('Problem displaying title montage for statistic_image object');
                    end
                end

                o2 = addblobs(o2, region(obj, 'noverbose'), 'nooutline', varargin{:});

            end

            drawnow

        end

        if ~do_multirow

            % Plot legend only for single-row, otherwise obscures some images

            o2 = legend(o2, varargin{:});  % pass through "noverbose" option

        end

        fig_handle = o2;
        clear o2

    case 'scnmontage'

        overlay = which('SPM8_colin27T1_seg.img');

        cl = cell(1, n);
        fig_handle = zeros(1, n);

        for i = 1:n

            % data from this image
            dat = image_obj.dat(:, i);

            % top and bottom 10%
            %dat(dat > prctile(dat, 10) & dat < prctile(dat, 90)) = 0;

            cl{i} = iimg_indx2clusters(dat, image_obj.volInfo);

            fig_handle(i) = montage_clusters(overlay, cl{i}, [2 2]);

            set(fig_handle, 'Name', sprintf('Montage %3.0f', i), 'Tag', sprintf('Montage %3.0f', i))

        end

    otherwise, warning(['Unknown input string option:' varargin{i}]);

end % switch

end % function

function [datvec, clim] = get_data_range(obj)
% see same function in render_on_surface, similar in render_blobs

clim = [];

if isa(obj, 'statistic_image')

    sig = logical(obj.sig);
    dat = obj.dat(sig);

else

    wh = obj.dat ~= 0 & ~isnan(obj.dat);
    dat = obj.dat(wh);

end

datvec = dat(:);

if any(isinf(datvec))
    warning('Some image values are Inf. Expect erratic behavior/errors.');
    whinf = isinf(datvec);
    datvec(whinf) = sign(datvec(whinf)) .* max(abs(datvec(~whinf)));
end

if isempty(clim)

    %clim = [min(datvec) max(datvec)];  % clim: data values for min and max, should become min/max colors

    if any(datvec < 0) && any(datvec > 0) % split colormap

        clim = double([prctile(datvec(datvec < 0), 10) prctile(datvec(datvec > 0), 90)]); % Match defaults for montage in render_blobs

    else
        % may need adjustment
        clim = [min(datvec) max(datvec)];

    end

    % if no variance, we have constant data values - special case.
    % reducing lower clim(1) will effectively map all data to max color value
    if abs(diff(clim)) < 100 * eps
        clim(1) = clim(2) - 1;
    end

end

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% Test Tx DMA data output
amplitude = 2^15; frequency = 20e6;
swv1 = dsp.SineWave(amplitude, frequency);
swv1.ComplexOutput = true;
swv1.SamplesPerFrame = 2^20;
swv1.SampleRate = 100e6;
y = swv1();

uri = 'ip:analog';
fc = 1e9;

%% Tx set up
tx = adi.AD9371.Tx('uri',uri);
tx.CenterFrequency = fc;
tx.EnableCustomProfile = true;
tx.CustomProfileFileName = 'profile_TxBW100_ORxBW100_RxBW100.txt';
tx.DataSource = 'DMA';
tx.EnableCyclicBuffers = true;
tx.AttenuationChannel0 = -10;
tx(y);

%% Rx set up
rx = adi.AD9371.Rx('uri',uri);
rx.CenterFrequency = fc;

%% Run
for k=1:20
    valid = false;
    while ~valid
        [out, valid] = rx();
    end
end
rx.release();
tx.release();

%% Plot
nSamp = length(out);
fs = tx.SamplingRate;
FFTRxData  = fftshift(10*log10(abs(fft(out))));
df = fs/nSamp;  freqRangeRx = (-fs/2:df:fs/2-df).'/1000;
plot(freqRangeRx, FFTRxData);
xlabel('Frequency (kHz)');ylabel('Amplitude (dB)');grid on;

author: analogdevicesinc


text:
function [out] = saturation_7(p1,p2,p3,p4,p5,S,In)
%saturation_7

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Saturation excess from a store with different degrees of
%               saturation (gamma function variant)
% Constraints:  f = 0, for x-p3 < 0
%               S >= 0      prevents numerical problems with integration
% @(Inputs):    p1   - scaling parameter [-]
%               p2   - gamma function parameter [-]
%               p3   - storage threshold for flow generation [mm]
%               p4   - absolute scaling parameter [mm]
%               p5   - linear scaling parameter [-]
%               S    - current storage [mm]
%               In   - incoming flux [mm/d]

out = integral(@(x)...
        1./(p1.*gamma(p2)).*(max(x-p3,0)./p1).^(p2-1).*exp(-1.*max(x-p3,0)./p1),...
        p5.*max(S,0)+p4,Inf).*In;

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function fighandle=newdfig(varargin)
% fighandle = newdfig(numfigs,figname)  create NEW Docked FIGures.
%   "fighandle" is a vector of handle numbers for the figures.
%   "numfigs" determines the number of figures to create.
%   "figname" is an optional string for each figure's name.
%   Examples:
%       newdfig('Title'), newdfig Title, newdfig
%       figh=newdfig(3); figh=newfig(2,'Title');

%   Copyright 2006 Mirtech, Inc.
%   created 08/27/2006  by Mirko Hrovat on Matlab Ver. 7.2
%   Mirtech, Inc.       email: mhrovat@email.com
numfigs=[];
ftitle=[];
if nargin~=0,
    for n=1:length(varargin),
        arg=varargin{n};
        if isnumeric(arg),
            numfigs=arg;
        elseif ischar(arg),
            ftitle=arg;
        end
    end  % for
end  % if nargin
if isempty(numfigs),    numfigs=1;  end
if isempty(ftitle),     ftitle='';  end

figh=zeros(1,numfigs);
for n=1:numfigs,
    figh(n)=figure('Name',ftitle,'WindowStyle','docked');
end

if nargout==1,
    fighandle=figh;
end


text:
function kernWriteToFID(kern, FID)

% KERNWRITETOFID Load from an FID written by the C++ implementation.
% FORMAT
% DESC loads in from a file stream the data format produced by
% C++ implementations.
% ARG kern : the kernel structure to write to the stream.
% ARG FID : the file ID from where the data is loaded.
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006, 2008
%
% SEEALSO : modelReadFromFID, kernCreate, kernReadParamsFromFID

% KERN

writeVersionToFID(FID, 0.2);
writeStringToFID(FID, 'baseType', 'kern');
writeStringToFID(FID, 'type', kern.type);
kernWriteParamsToFID(kern, FID);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function V = LP10cvx(K, P, model, epsilon)
% CPLEX implementation of LP-9 for input sets K, P (see FASTCORE paper)
%
% USAGE:
%
%    V = LP9cvx(K, P, model, epsilon)
%
% .. Authors: - Nikos Vlassis, Maria Pires Pacheco, Thomas Sauter, 2013 LCSB / LSRU, University of Luxembourg

scalingfactor = 1e3;

V = [];
if isempty(P) || isempty(K)
    return;
end

np = numel(P);
nk = numel(K);
n = size(model.S,2);

cvx_begin

  variable v(n);
  variable z(np);

  minimize( ones(1,np) * z );

  z>=0;
  v(P)>=-z;
  v(P)<=z;

  v(K)>=epsilon*scalingfactor;

  model.S*v==0;

  v>=model.lb*scalingfactor;
  v<=model.ub*scalingfactor;

cvx_end

V = v;

author: opencobra

repo: cobratoolbox


text:
function [grad_W_rnn, grad_W_emb, emb_indices] = srcCharLayerBackprop(W_rnn, charData, charGrad)
% Backprop for char layer from word gradients to chars.
% Input:
%   W_rnn: recurrent connections of multiple layers, e.g., W_rnn{ll}.
%
% Thang Luong @ 2015, <lmthang@stanford.edu>

  if charData.numRareWords == 0
    grad_W_rnn = [];
    grad_W_emb = [];
    emb_indices = [];
    return;
  end

  assert(length(charGrad.indices) == charData.numRareWords);

  params = charData.params;
  if params.assert
    assert(isequal(sort(charData.rareWordMap(charGrad.indices))', 1:charData.numRareWords));
  end

  topGrads = cell(charData.maxLen, 1);
  topGrads{end} = charGrad.embs(:, charData.rareWordMap(charGrad.indices));

  % init state
  zeroBatch = zeroMatrix([params.lstmSize, params.curBatchSize], params.isGPU, params.dataType);
  zeroState = cell(params.numLayers, 1);
  zeroGrad = cell(params.numLayers, 1);
  for ll=1:params.numLayers % layer
    zeroState{ll}.h_t = zeroBatch;
    zeroState{ll}.c_t = zeroBatch;
    zeroGrad{ll} = zeroBatch;
  end

  [~, ~, grad_W_rnn, grad_W_emb, emb_indices, ~, ~] = rnnLayerBackprop(W_rnn, charData.states, zeroState, ...
  topGrads, zeroGrad, zeroGrad, charData.batch, charData.mask, charData.params, charData.rnnFlags, [], [], [], []);
end

author: lmthang


text:
function sz = tsize(a,idx)
%TSIZE Tensor size of sptenmat.
%
%   D = TSIZE(X) returns the size of the tensor being stored as a
%   matrix.
%
%   M = TSIZE(X,DIM) returns the length of the dimension(s) specified
%   by DIM.  For example, SIZE(X,1) returns the size of the first
%   dimension of the tensor.
%
%   See also SPTENMAT, SPTENMAT/SIZE.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if isempty(a.tsize)
    sz = [];
    return;
end

if exist('idx', 'var')
    sz = a.tsize(idx);
else
    sz = a.tsize;
end

return;

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function str= str_vec2str(v, fmt, inbetween)
%str= str_vec2str(v, <format, inbetween=', '>)
%
% v can be a vector of real numbers or a cell array of strings.
% the default for format is '%d' in the former and '%s' in the latter case.

%narginchk(1,3)
nargchk(1, 3, nargin);
misc_checkType(v,'DOUBLE[-]|CELL{CHAR}');
misc_checkTypeIfExists('fmt','CHAR');
misc_checkTypeIfExists('inbetween','CHAR');

if ~exist('fmt', 'var') || isempty(fmt),
  if iscell(v) || ischar(v),
    fmt='%s';
  else
    fmt='%d';
  end
end
if ~exist('inbetween', 'var'), inbetween=', '; end

if isempty(v),
  str= '[]';
elseif iscell(v),
  isprint= 1:length(v);
  for vi= 1:length(v),
    if ischar(v{vi}),
      fstr{vi}= ['%s' inbetween];
    elseif ~iscell(v{vi}) && ~isstruct(v{vi}),
      fstr{vi}= repmat([fmt inbetween], 1, length(v{vi}));
    else
      isprint= setdiff(isprint, vi,'legacy');
    end
  end
  if isempty(isprint), str=''; return; end
  v= {v{isprint}};
  fstr= {fstr{isprint}};
  if length(isprint)>1,
    str= [sprintf([fstr{1:end-1}], v{1:end-1}) ...
          sprintf(fstr{end}(1:end-length(inbetween)), v{end})];
  else
    str= sprintf(fstr{end}(1:end-length(inbetween)), v{end});
  end
else
  if length(v)==1 || ischar(v),
    str= sprintf(fmt, v);
  else
    str= [sprintf([fmt inbetween], v(1:end-1)) ...
          sprintf(fmt, v(end))];
  end
end

author: PatternRecognition

repo: OpenBMI


text:
function [isLFP] = bz_isLFP(lfp)
% USAGE
% [] = bz_isLFP(lfp)
%
% INPUT
%       lfp   - struct with the following fields
%                   .lfp
%                   .timestamps
%                   .samplingRate
%                   .channels
%
%
% OUTPUT
%      logical true if struct meets lfp criteria, false if otherwise
%
% written by david tingley, 2017

if isfield(lfp,'data') && isfield(lfp,'timestamps') && isfield(lfp,'samplingRate') && isfield(lfp,'channels') % check that fields exist
     if ismatrix(lfp.data) && isvector(lfp.timestamps) && isnumeric(lfp.samplingRate) && isvector(lfp.channels)
         isLFP = true;
     else
         warning('one of the required fields for a lfp type is not formatted correctly')
         isLFP = false;
     end
else
    warning('one of the required fields for a lfp type does not exist')
    isLFP = false;
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ h, hmin, hmax ] = p20_stepsize ( option )

%*****************************************************************************80
%
%% P20_STEPSIZE returns step sizes for problem 20.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Output, real H, HMIN, HMAX, suggested values for the
%    initial step, the minimum step, and the maximum step.
%
  h =    0.00250;
  hmin = 0.01000;
  hmax = 0.08000;

  return
end


text:
% THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION IS RELEASED "AS IS."  THE U.S. GOVERNMENT MAKES NO WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, CONCERNING THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION, INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT WILL THE U.S. GOVERNMENT BE LIABLE FOR ANY DAMAGES, INCLUDING ANY LOST PROFITS, LOST SAVINGS OR OTHER INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE, OR INABILITY TO USE, THIS SOFTWARE OR ANY ACCOMPANYING DOCUMENTATION, EVEN IF INFORMED IN ADVANCE OF THE POSSIBILITY OF SUCH DAMAGES.
%
% change log:
%   040204 smw stub file
%
% file: get_r_ci_tl.m
%
function varargout = get_r_ci_tl(varargin)
  error('get_r_ci_tl is not part of current release.');
return;


text:
function [ y, next ] = p03_equil ( neqn, next )

%*****************************************************************************80
%
%% P03_EQUIL returns equilibrium solutions of problem p03.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, integer NEXT, the index of the previous
%    equilibrium, which should be 0 on first call.
%
%    Output, real Y(NEQN), the "next" equilibrium solution, if any.
%
%    Output, integer NEXT, the index of the current equilibrium,
%    or 0 if there are no more.
%
  if ( next == 0 )
    next = 1;
    y(1:neqn,1) = 0.0;
  else
    next = 0;
    y = [];
  end

  return
end


text:
% Calculate border score for a firing map
%
% Calculates a border score for a firing rate map according to the article "Representation of Geometric
% Borders in the Entorhinal Cortex" by Solstad et. al. (Science 2008).
% Border score ranges from -1 to +1 with +1 being "a perfect border cell". If the firing map contains
% no firing fields, then the returned score is -1.
% The score reflects not only how close a field is to a border and how big the coverage of this field is,
% but also it reflects spreadness of a field. The best border score (+1) will be calculated for a thin
% line (1 px, bin) that lies along the wall and fully covers this wall.
%
%  USAGE
%   score = analyses.borderScore(map, fieldsMap, fields, <options>)
%   map         A 2D firing map, not a struct you get from analyses.map, but actual values.
%   fieldsMap   Matrix of the same size as map with only non-zero elements
%               belonging to detected fields. This is an output from
%               function analyses.placefield.
%   fields      Array of structures with information about detected fields.
%               This is an output from function analyses.placefield.
%   <options>   Optional list of property-value pairs (see table below)
%
%   ==============================================================================================
%    Properties    Values
%   ----------------------------------------------------------------------------------------------
%    'searchWidth'  If map is not perfect, but contains NaN values along borders, then
%                   search for border pixels can have NaNs. To mitigate this, we check
%                   searchWidth rows/columns near border and if the closest to the border pixel
%                   equals to NaN, we search for first non-NaN value in searchWidth rows-columns.
%                   This argument is optional and default value is 8.
%    'walls'        Definition of walls along which the border score is calculated. Provided by
%                   a string which contains characters that stand for walls:
%                   T - top wall (we assume the bird-eye view on the arena)
%                   R - right wall
%                   B - bottom wall
%                   L - left wall
%                   Characters are case insensitive. Default value is 'TRBL' meaning that border
%                   score is calculated along all walls. Any combination is possible, e.g.
%                   'R' to calculate along right wall, 'BL' to calculate along two walls, e.t.c.
%   ==============================================================================================
%   score       Border score. Ranges from -1 to +1 (if there are no fields,
%               then score is -1).
%
%  SEE
%   See analyses.placefield
%
function score = borderScore(map, fieldsMap, fields, varargin)
    score = -1;

    if ~isempty(fields)
        coverage = analyses.borderCoverage(fields, varargin{:});

        fieldsFiringMap = fieldsMap;
        fieldsFiringMap(fieldsFiringMap > 0) = 1; % different fields are assigned different numbers,
                                                  % make them all 1.
        fieldsFiringMap = map .* fieldsFiringMap;

        fdist = weighted_firing_distance(fieldsFiringMap);
        score = (coverage - fdist)/(coverage + fdist);
    end
end

function wfd = weighted_firing_distance(map)
    % normalized firing map. Normalization is done by
    % the sum of firing rates of all pixels belonging to all fields
    map = map/nansum(map(:));

    [ly, lx] = size(map);
    [mx, my] = meshgrid(1:lx, 1:ly);

    % alternative code for distance matrix is:
    % map = zeros(ly, lx);
    % map(1, :) = 1;
    % map(end, :) = 1;
    % map(:, 1) = 1;
    % map(:, end) = 1;
    % D = bwdist(map);
    % but this code is slower.

    distance_matrix = min(min(my,mx), min(flipud(my), fliplr(mx)));

    wfd = nansum(nansum(map .* distance_matrix));

    % normalization by half of the smallest arena size min(ly, lx)/2.
    wfd = (2 * wfd) / min(ly, lx);
end

author: kavli-ntnu


text:
%Author : Athi Narayanan S
%M.E, Embedded Systems,
%K.S.R College of Engineering
%Erode, Tamil Nadu, India.
%http://sites.google.com/site/athisnarayanan/

function mergedBlks = merge_blocks(data)

s = size(data);
mergedBlks = zeros((s(3) * s(1)), (s(4) * s(2)));

%Merging the Blocks
for i = 1 : s(3)
    for j = 1 : s(4)
        localBlk = data(:,:,i,j);
        rowStart = ((i - 1) * s(1)) + 1;
        rowEnd = i * s(1);
        colStart = ((j - 1) * s(2)) + 1;
        colEnd = j * s(2);
        mergedBlks(rowStart:rowEnd, colStart:colEnd) = localBlk;
    end
end


text:
function gmm = mixtureOfGaussian(data, K, niter, fullCov, varargin)
% Mixture of K gaussians
% data(ndata, nvariables) - the data
% K - the number of gaussians in the mixture
% niter - the number of iterations
% fullCov - 1 if a full covariance matrix should be used, otherwise
%                   a diagonal matrix will be used
% varargin{1} - weights for the data (optional)
% gmm.mu{K} - the means of the mixtures
% gmm.sigma{K} - the variance matrices
% gmm.priors(K) - the mixture priors

[ndata, ndim] = size(data);

if ndata == 0
    disp('warning no data')
    gmm.mu = {};
    gmm.sigma = {};
    gmm.priors = [];
    return;
end

if length(varargin)==0
    w = ones(ndata, 1);
else
    w = varargin{1}(:);
end

w = w / sum(w) * ndata;

cumw = cumsum(w / sum(w));

% initialize mu by randomly selected data values
r = rand(K, 1);
for k = 1:K
    for i = 1:ndata
        if cumw(i)>=r(k)
            mu{k} = data(i, :);
            break;
        end
    end
end

% initialize all covariances to be covariance of data
%if fullCov
%    sigma(1:K) = {cov(data)};
%else
%    sigma(1:K) = {diag(var(data))};
%end
sigma(1:K) = {diag(ones(1, ndim))};

priors = ones(1, K) / K;

for iter = 1:niter

    % get posteriors: P(k | data) = P(data|k)P(k) / sum_k(P(data|k))
    pC = zeros(ndata, K);
    for k = 1:K
        invSigma = inv(sigma{k});
        for i = 1:ndata
            pC(i, k) = 1 / (sqrt(2*pi)*sqrt(det(sigma{k}))) * ...
                exp(-1/2 * (data(i, :)-mu{k}) * invSigma * (data(i, :)-mu{k})') ...
                * priors(k);
        end
    end

    if 0
    for k = 1:K
        disp(['k: ' num2str(k)])
        disp(['mu: ' num2str(mu{k})])
        disp(['sigma: ' num2str(det(sigma{k}))])
        disp(['prior: ' num2str(priors(k))])
        disp(['min: ' num2str(min(pC(:, k)))])
        disp(['max: ' num2str(max(pC(:, k)))])
    end
    end

    %disp(num2str(priors))

    pC = pC + eps; % avoid divide by zero error

    sumPc = sum(pC, 2);
    for k = 1:K
        pC(:, k) = pC(:, k) ./ sumPc .* w;
    end
    kweight = sum(pC, 1);

    priors = kweight / sum(kweight);

    % mean estimation
    for k = 1:K
        sumkX = pC(:, k)' * data;
        mu{k} = sumkX / kweight(k);
    end

    % covariance matrix estimation
    if fullCov
        for k = 1:K
            diffs = data - (ones(ndata, 1) * mu{k});
            diffs = diffs .* (sqrt(pC(:,k))*ones(1, ndim));
            sigma{k} = (diffs' * diffs) / kweight(k);
        end
    else
        for k = 1:K
            diffs = data - (ones(ndata, 1) * mu{k});
            sigma{k} = diag(sum((diffs.*diffs).*(pC(:, k)*ones(1, ndim)), 1) ./ kweight(k));
        end
    end
end

gmm.mu = mu;
gmm.sigma = sigma;
gmm.priors = priors;

author: zouchuhang

repo: LayoutNet


text:
function x=v_rsfft(y,n)
%V_RSFFT    fft of a real symmetric spectrum X=(Y,N)
% Y is the "first half" of a symmetric real input signal and X is the
% "first half" of the symmetric real fourier transform.
% If the length, N, of the full signal is even, then the "first half"
% contains 1+N/2 elements (the first and last are excluded from the reflection).
% If N is odd, the "first half" conatins 0.5+N/2 elements and only the first
% is excluded from the reflection.
% If N is specified explicitly, then Y will be truncated of zero-padded accordingly.
% If N is omitted it will be taken to be 2*(length(Y)-1) and is always even.
%
% If Y is a matrix, the transform is performed along each column
%
% The inverse function is y=v_rsfft(x,n)/n

% Could be made faster for even n by using symmetry

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: v_rsfft.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 2 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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~isreal(y) error('RSFFT: Input must be real'); end
fl=size(y,1)==1;
if fl y=y(:); end
[m,k]=size(y);
if nargin<2 n=2*m-2;
else
    mm=1+fix(n/2);
    if mm>m y=[y; zeros(mm-m,k)];
    elseif mm<m y(mm+1:m,:)=[];
    end
    m=mm;
end
x=real(fft([y;y(n-m+1:-1:2,:)]));
x(m+1:end,:)=[];
if fl x=x.'; end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function[varargout]=blurspec(varargin)
%BLURSPEC  Returns the blurred and aliased spectrum given the autocovariance.
%
%   BLURSPEC is used to rapidly compute blurred, aliased, and other
%   modified versions of a spectrum from a known autocovariance.
%
%   Performing these calculations in the time domain, making use of
%   Fourier relationships between modified versions of the autocovariance
%   and spectrum, is much faster than working in the frequency domain.
%   __________________________________________________________________
%
%   Blurred and aliased spectra
%
%   [F,S]=BLURSPEC(R) inverse Fourier transforms a given *one-sided*
%   autocovariance R to obtain a spectrum S that incorporates the effects
%   of aliasing, as well as blurring by the default taper, the 'boxcar'.
%
%   Here F is an array of frequencies and S is the one-sided Fourier
%   spectrum.  If R is length N, the output arrays will be length (N/2+1)
%   if N is even and (N+1)/2 if N is odd.
%
%   [F,SPP,SNN]=BLURSPEC(R) for complex-valued R returns SPP and SNN, the
%   positive and negative rotary spectra, respectively.
%
%   BLURSPEC(DT,R) optionally uses the sample interval DT in computing the
%   frequency array F, and in setting the spectral levels.
%   __________________________________________________________________
%
%   Tapered and aliased spectra
%
%   BLURSPEC(R,'tapered',TAPER) incorporates the spectral smoothing due to
%   the use of data TAPER, rather than the boxcar taper, in addition to the
%   effects of aliasing.  TAPER must be the same length as R.
%
%   Note that BLURSPEC(R,'tapered',[]) simply returns the blurred spectrum,
%   as an empty taper is taken to indicate the periodogram.
%
%   BLURSPEC(R,'window',WIN) uses the pre-computed window that is to
%   multiply the autocovariance function WIN.  This is the half of the
%   sequence obtained by convolving the taper with itself.  This version
%   is primarily used for speed in an internal call from MATERNFIT.
%   __________________________________________________________________
%
%   Aliased-only spectrum
%
%   BLURSPEC(R,'aliased') computes an approximation to *aliased* spectrum,
%   without blurring.  This approximation will be accurate to the extent
%   that R has decayed to zero by the end of its duration.
%
%   This is much faster than explicitly summing over aliased frequencies.
%   __________________________________________________________________
%
%   Differenced spectra
%
%   BLURSPEC(R,'difference') returns the blurred and aliased spectrum of
%   the first forward difference of the process with autocovariance R.  In
%   this case the output will be length N-1.
%
%   BLURSPEC(R,'seconddifference') returns the blurred and aliases spectrum
%   of the second forward difference of the process with autocovariance R.
%   The output will be length N-2.
%
%   These options can be combined with the 'taper' and 'aliased' options
%   described above.
%   __________________________________________________________________
%
%   'blurspec --t' runs some tests.
%
%   Usage: [f,S]=blurspec(R);
%          [f,Spp,Snn]=blurspec(dt,R);
%          [f,Spp,Snn]=blurspec(dt,R,'tapered',TAPER);
%          [f,Spp,Snn]=blurspec(dt,R,'aliased');
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2015--2020 J.M. Lilly and A.M. Sykulski
%                                 --- type 'help jlab_license' for details

if strcmp(varargin{1}, '--t')
    blurspec_test,return
end

dt=1;
if length(varargin{1})==1
    dt=varargin{1};
    varargin=varargin(2:end);
end
R=varargin{1};

ver='standard';
str='blurred';
psi=[];

varargin=varargin(2:end);
for i=1:2
    if length(varargin)>0
        if ischar(varargin{end})
            if strcmpi(varargin{end}(1:3),'sec')||strcmpi(varargin{end}(1:3),'dif')||strcmpi(varargin{end}(1:3),'sta')
                ver=varargin{end};
            elseif strcmpi(varargin{end}(1:3),'ali')
                str=varargin{end};
            end
            varargin=varargin(1:end-1);
        end
    end
    if length(varargin)>1
        %length(varargin)
        if ischar(varargin{end-1})
            str=lower(varargin{end-1});
            psi=varargin{end};
            varargin=varargin(1:end-2);
        end
    end
end

if isempty(psi)&&strcmpi(str(1:3),'tap')
    str='blurred';
end

if strcmpi(ver(1:3),'dif')
    R=frac(1,dt)*(2*R(1:end-1,:)-R(2:end,:)-[conj(R(2,:));R(1:end-2,:)]);
elseif strcmpi(ver(1:3),'sec')
    R=frac(1,dt)*(2*R(1:end-1,:)-R(2:end,:)-[conj(R(2,:));R(1:end-2,:)]);
    R=frac(1,dt)*(2*R(1:end-1,:)-R(2:end,:)-[conj(R(2,:));R(1:end-2,:)]);
end

N=size(R,1);
if strcmpi(str(1:3),'win')
    R=R.*vrep(psi,size(R,2),2);
elseif strcmpi(str(1:3),'tap')
    if length(psi)~=size(R,1)
        error('Sizes of R and taper PSI do not match.')
    end
    win=conv(psi,psi);
    win=win(end-N+1:end);
    R=R.*vrep(win,size(R,2),2);
elseif strcmpi(str(1:3),'blu')
    tri=[N:-1:1]'./N;
    if size(R,2)==1
        R=R.*tri;
    else
        R=R.*vrep(tri,size(R,2),2);
    end
end  %If str='aliased', do nothing
%figure,plot(tri)

R(1,:)=R(1,:)./2;  %Don't forget to divide first element by two
S=dt*2*real(fft(R));
S=abs(S);  %Sometimes there are small negative parts after blurring

%Note, this is always correct for both even and odd length time series
omega=fourier(N);
Spp=S(1:length(omega),:);
Snn=[S(1,:);S(end:-1:end-length(omega)+2,:)];
%Snn=flipud([S(end-length(omega)+2:end,:);S(1,:)]);  %Same but slower

varargout{1}=omega./dt;
varargout{2}=Spp;
varargout{3}=Snn;

function[]=blurspec_test

N=1000;
alpha=1.5;
h=1;
A=1;

[tau,R]=materncov(1,N,A,alpha,h);
[f,Spp,Snn]=maternspec(1,N,A,alpha,h);
tic;[f,Spp2,Snn2]=blurspec(R,'aliased');etime1=toc;

maternspec_spec=@(omega,A,omegao,H,alpha)(frac(H.^(2*alpha-1),materncfun(alpha))*frac(A.^2,((omega-omegao).^2+H.^2).^alpha));
tic;
M=100;
[Sppa,Snna]=vzeros(length(fourier(N)),2*M+1);
for m=-M:M
    Sppa(:,m+M+1)=maternspec_spec(fourier(N)+2*pi*m,A,0,h,alpha);
end
Snna=Sppa;
etime2=toc;

vsum(Sppa,Snna,2);
%figure,plot(f,[Spp2 Spp2],'k','linewidth',2),hold on,plot(f,[Sppa Snna],'r')
%figure,plot(f,Spp2-Sppa,'k','linewidth',2),hold,plot(f,Snn2-Snna,'r','linewidth',2)

b1=aresame(Spp2,Sppa,2e-6);
b2=aresame(Snn2,Snna,2e-6);

reporttest('BLURSPEC aliasing with covariance inversion matches direct calculation to 2e-6, even N',b1&&b2)
disp(['BLURSPEC aliasing with covariance inversion was ' num2str(etime2/etime1) ' times faster than direct calculation.'])

N=999;
omegao=1/4;
[tau,R]=materncov(1,N,A,alpha,h,omegao);
[f,Spp,Snn]=maternspec(1,N,A,alpha,h,omegao);
tic;[f,Spp2,Snn2]=blurspec(R,'aliased');etime1=toc;

tic;
M=200;
[Sppa,Snna]=vzeros(length(fourier(N)),2*M+1);
for m=-M:M
    Sppa(:,m+M+1)=maternspec_spec(fourier(N)+2*pi*m,A,omegao,h,alpha);
    Snna(:,m+M+1)=maternspec_spec(fourier(N)+2*pi*m,A,-omegao,h,alpha);
end
etime2=toc;

vsum(Sppa,Snna,2);
b1=aresame(Spp2,Sppa,1e-6);
b2=aresame(Snn2,Snna,1e-6);

reporttest('BLURSPEC aliasing with covariance inversion matches direct calculation to 1e-6, odd N, frequency shift',b1&&b2)
disp(['BLURSPEC aliasing with covariance inversion was ' num2str(etime2/etime1) ' times faster than direct calculation.'])

x=[10 1.1 0.1];
[tau,acv]=materncov(1,N,x(1),x(2),x(3)); % autocovariance sequence
S=abs(real(2*fft(acv.*(1-([0:N-1]')/N))-acv(1))); % blurred spectrum

Spp2=S(1:floor(N/2)+1);
Snn2=[S(1);S(end:-1:ceil(N/2)+1)];
[tau,R]=materncov(1,N,x(1),x(2),x(3));
[f,Spp,Snn]=blurspec(R);

%[f,Spp,Snn]=maternspec(N,x(1),x(2),x(3),'blurred');

reporttest('BLURSPEC blurring matches alternate version for Matern', aresame(Spp,Spp2,1e-8)&&aresame(Snn,Snn2,1e-8))

% N=1000;
% alpha=10;
% h=1/10;
% sigma=7;
% z=maternoise(1,[N 10000],sigma,-1/2,h,0,alpha);
% psi=sleptap(size(z,1),3,1);
% [f,spp,snn]=mspec(1,z,psi);
%
%
% [tau,R]=materncov(1,N,sigma,-1/2,h,0,alpha);
% [f,Spp,Snn]=blurspec(R,'taper',psi);
% figure,plot(f,vmean([spp snn],2)),hold on,plot(f,Spp)
%
% N=1001;
% z=maternoise(1,[N 10000],sigma,-1/2,h,0,alpha);;
% psi=sleptap(size(z,1),3,1);
% [f,spp,snn]=mspec(1,z,psi);
%
%
% [tau,R]=materncov(1,N,sigma,-1/2,h,0,alpha);
% [f,Spp,Snn]=blurspec(R,'taper',psi);
% %figure,plot(f,vmean([spp snn],2)),hold on,plot(f,Spp)
%
%
% N=1000;
% alpha=1.5;
% h=1;
% A=1;
% psi=sleptap(size(z,1),3,1);
%
% [f,Spp,Snn]=maternspec(dt,N,A,alpha,h/dt);
% [tau,R]=materncov(dt,N,A,alpha,h);
% tic;[f,Spp2,Snn2]=blurspec(dt,R,'tapered',psi);
% tic;[f,Spp3,Snn3]=blurspec(dt,R);
% figure,plot(f,[Spp Spp2 Spp3])
%
%
% dt=7;
% z=maternoise(dt,[N 1000],A,alpha,h/dt);
% psi=sleptap(size(z,1),3,1);
% [f,spphat,snnhat,spn]=mspec(dt,z,psi);
% vmean(spphat,snnhat,2);
%
%

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function allowaxestogrow(f)
%
% Call allowaxestogrow(f) on a figure f with multiple sub-axes.
%
% Then, whenever one of the axes is clicked, it zooms to
% the size of the window (or a preset user-defined size).
% When clicked again, the axis returns to its original position.
%
% Useful when you have a number of tiny subplots in a figure window,
% and would like to be able to visually inspect them without squinting.
%
% To set the zoomed (large) size for an axis ax:
%   setappdata(ax,'LastPosition',[LEFT BOTTOM WIDTH HEIGHT]);
% For colorbars, [0 0 0.85 1] usually works well.
%
% Version history:
%
% 0.40: Use get/setappdata instead of userdata. Use only one button-down
%       function. Add some extra links between axes for extensibility.
%       Suggestions by Jonas. 2009.05
% 0.32: Ensure invisible patch is created with correct figure parent. 2009.05
% 0.31: Try to preserve plot order better; reorder plots upon shrinking
%       so that overlay axes, etc. are better respected. Still not perfect --
%       ideally we'd preserve the original exact plot order. 2009.05
%       Preserve plot aspect ratio when creating invisible axes, following
%       suggestion by Jonas. 2009.04
%
% 0.30: Much simplified by using invisible patches to trap clicks.
%       Now requires Matlab's OpenGL renderer which may cause side effects
%       on some systems.
%       Also, allows user to set predefined size for zoomed plot.
%       Among other things, this makes colorbars finally work well. 2009.02
%
% 0.22: add workaround for reported bug when axis type is a cell array.
%       also add click support for patches.  2007.06
% 0.21: add some support for colorbars 2007.06.
% 0.20: modified to support clicks on images and lines following suggestion
%       from Jessica Louie 2007.05
% 0.11: Added code from Iram Weinstein to keep legends
%       from disappearing. 2005.08
% 0.10: Written by Matthew Caywood, <caywood@phy.ucsf.edu>, 2005.06
%
% Known bugs:
%
% Possible error if axes are plotted on top of panels.
% I don't use panels, so I'm not currently supporting them.
% Bugfixes are welcome, though. (reported by Scott Hirsch)
%
% With multiple axes in a subplot, clicks will sometimes magnify the wrong
% one, making the plot order wonky. For now, just keep clicking and it
% should fix itself.

if nargin==0, f=gcf; end

set(f,'Renderer','OpenGL');

% MSC necessary to avoid an OpenGL renderer bug where axes are flipped
% If you're having rendering problems, comment out the following line:
opengl('OpenGLBitmapZbufferBug',0)

% Cover all subplot axis children of figure f with an identically sized axis
% containing an invisible patch which triggers an event when clicked.
% Do not add to legends or colorbars.
ch = get(f,'Children');
for i = 1:length(ch)
    ax = ch(i);
    if (strcmp(get(ax,'Type'),'axes')) && ...
       ~strcmp(get(ax,'Tag'),'legend') && ...
       ~strcmp(get(ax,'Tag'),'Colorbar')

        % invisible invisax contains patch p which traps clicks,
        invisax = axes('Position',get(ax,'Position'), ...
                       'PlotBoxAspectRatioMode',get(ax,'PlotBoxAspectRatioMode'),...
                       'DataAspectRatioMode',get(ax,'DataAspectRatioMode'),...
                       'NextPlot','add');
        % invisax contains a link to ax in appdata
        setappdata(invisax,'VisibleAxes',ax);
        setappdata(invisax,'State',0); % 1 = enlarged, 0 = normal/shrunk
        axis off;

        p = patch([0 0 1 1],[0 1 1 0],[0 0 0 0],'w');
        set(p,'FaceAlpha',0);
        set(p,'EdgeAlpha',0);
        set(p,'ButtonDownFcn',@togglesubplot);

        % tag the invisible axes so that they can be manipulated as necessary
        set(p,'Tag','catchPatch')
        setappdata(ax,'catchPatchHandle',p);
        setappdata(ax,'InvisibleAxes',invisax);

        % ensure invisible patch gets put on right figure
        % MSC should not be necessary -- may work around
        % bug when creating figures quickly, at least in R14
        set(invisax,'Parent',f);
    end
end

% -------------

function togglesubplot(src,eventdata) %#ok

% default size of enlarged axes, modify as necessary
presetaxes = [0.05 0.05 0.9 0.9];

% patch p was clicked, get invisax and ax
p = src;
invisax = get(p,'Parent');
ax = getappdata(invisax,'VisibleAxes');
state = getappdata(invisax,'State');

% get current & last position
currentpos = get(ax,'Position');
futurepos = getappdata(ax,'LastPosition');
if isempty(futurepos), futurepos = presetaxes; end

% swap current & last position
set(ax,'Position',futurepos);
set(invisax,'Position',futurepos);

setappdata(ax,'LastPosition',currentpos);

% get order of parent figure's children
f = get(ax,'Parent');
ch = get(f,'Children');

if state == 1
    % shrink
    % send invisax and ax to back
    newch = [ch(ax ~= ch & invisax ~= ch) ; invisax ; ax];
    set(f,'Children',newch);
else
    % enlarge
    % change order of parent figure's children so clicked axis
    % is drawn on top.  But keep legend uppermost, if present.

    % Find a legend for the current subplot, if it exists
    ax = findobj(ax,'type','axes');
    leg = find_legend(ax);

    % move invisax and ax to top, then move legend if present
    newch = [invisax ; ax ; ch(ax ~= ch & invisax ~= ch)];
    if ~isempty(leg)
        newch = [leg ; newch(leg ~= newch)];
    end
    set(f,'Children',newch);
end

setappdata(invisax,'State',~state); % toggle state

%----------------------------------------------------%

function leg = find_legend(ha)
% returns handle to legend
% from legend.m

fig = ancestor(ha,'figure');
ax = findobj(fig,'type','axes');
leg=[];
k=1;
while k<=length(ax) && isempty(leg)
    if islegend(ax(k))
        hax = handle(ax(k));
        if isequal(double(hax.axes),ha)
            leg=ax(k);
        end
    end
    k=k+1;
end

%----------------------------------------------------%

function tf=islegend(ax)
% from legend.m

if length(ax) ~= 1 || ~ishandle(ax)
    tf=false;
else
    tf=isa(handle(ax),'scribe.legend');
end

author: Sable

repo: mcbench-benchmarks


text:
function [lx, ly, lz, varargin] = parseGridArgs3d(varargin)
%PARSEGRIDARGS3D  Extract or compute position vectors for meshgrid
%
%   [LX, LY] = parseGridArgs3d(LX, LY)
%   simply returns LX and LY.
%
%   [LX, LY] = parseGridArgs3d([NX NY])
%   assumes LX = 1:NX and LY = 1:NY
%
%   [LX, LY] = parseGridArgs3d([NX NY])
%   assumes LX = 1:NX and LY = 1:NY
%
%
%   Example
%   [lx, ly, lz] = parseGridArgs3d([100 100 100]);
%   [lx, ly, lz] = parseGridArgs3d(1:2:50, 1:4:100, 1:25);
%   [lx, ly, lz, varargin] = parseGridArgs3d(varargin{:});
%   for usage within another function.
%
%   See also
%     parseGridArgs

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2009-05-29,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

% If empty arguments, return default values
if isempty(varargin)
    lx = 1:100;
    ly = 1:100;
    lz = 1:100;
    return
end

var = varargin{1};

% case of a 2x3 matrix with starting position, increment, end position for
% each coordinate
if all(size(var) > [2 2])
    lx = var(1,1):var(1,2):var(1,3);
    ly = var(2,1):var(2,2):var(2,3);
    lz = var(3,1):var(3,2):var(3,3);
    varargin(1) = [];
    return;
end

% first argument contains maximal position for each coordinate
if all(size(var) == [1 3])
    lx = 1:var(1);
    ly = 1:var(2);
    lz = 1:var(3);
    varargin(1) = [];
    return;
end

% first and second arguments contain vector for each coordinate
% respectively
if length(varargin) > 2
    lx = varargin{1};
    ly = varargin{2};
    lz = varargin{3};
    varargin(1:3) = [];
    return
end

% otherwise, throws error
error('Error in parsing grid arguments');

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
%% t_meshShow
%
% This script loads a mesh and calls mrMesh to show it.
%
%
% See also:  t_meshFromClass
%
% (c) VISTASOFT, Stanford

%% Load an existing mesh and show it
tmp = load(fullfile(mrvDataRootPath,'anatomy','anatomyNIFTI','leftMesh.mat'));
msh = tmp.msh; clear tmp
msh = meshSet(msh,'actor',33);

windowID = 1000;
msh = meshSet(msh,'window id',1000);
meshVisualize(msh);

%% End

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function linplus_test25 ( )

%*****************************************************************************80
%
%% TEST25 tests R8GB_PRINT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  m = 8;
  n = 10;
  ml = 1;
  mu = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST25\n' );
  fprintf ( 1, '  For a general banded matrix,\n' );
  fprintf ( 1, '  R8GB_PRINT prints the matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix rows M =      %d\n', m );
  fprintf ( 1, '  Matrix columns N =   %d\n', n );
  fprintf ( 1, '  Lower bandwidth ML = %d\n', ml );
  fprintf ( 1, '  Upper bandwidth MU = %d\n', mu );

  a = r8gb_indicator ( m, n, ml, mu );

  r8gb_print ( m, n, ml, mu, a, '  The R8GB matrix:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function object=PASemptyobject
  object.label='';
  object.orglabel='';
  object.bbox=[];
  object.polygon=[];
  object.mask='';
  object.class='';
  object.view='';
  object.truncated=false;
  object.difficult=false;
return

author: quantombone

repo: exemplarsvm


text:
function r = emptyfh()
r = @(varargin)[];

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function writePOLY_tetgen(filename,V,F,H,varargin)
  % WRITEPOLY_TETGEN prints vertices and planar facets to a .poly file for
  % tetgen
  %
  % writePOLY_tetgen(filename,V,F,H)
  %
  % Input
  %   filename:  name of output file as string (caution! will clobber
  %                    existing)
  %   V  #V by dim list of vertex positions
  %   F  #F struct containing polygon information arrays
  %     .facets  a #facets list of facets, each facet is a again a list of
  %       polygons
  %      ** Note: contrary to writePOLY_pyramid, here facets index V directly
  %     .boundary_markers a #facets list of boundary_markers
  %     .holes  a #facets list of holes, each holes entry is again a list for
  %       each facet
  %    or
  %   F  #F by uniform_facet_size list of facets into V
  %   H  #H by dim list of volume hole positions
  %   Optional:
  %     'BoundaryMarkers' followed by #facets list of boundary markers
  %     'BoundaryMarkersNodes' followed by #facets list of boundary markers
  %     'RegionList' followed by #regions list of region vertex positions and (optional) region numbers and attributes
  %
  % Example:
  %   % constrained Delaunay tetrahedralization of a triangle mesh
  %   [V,F] = load_mesh('gargoyle.obj');
  %   % constrained Delaunay tetrahedralization of just points
  %   D = DelaunayTri(V);
  %   % Faces on boundary of convex hull
  %   BF = boundary_faces(D.Triangulation);
  %   % avoid dupicate faces
  %   FmBF = setdiff(sort(F,2),sort(BF,2),'rows');
  %   Facets = [];
  %   Facets.facets = mat2cell([BF;FmBF],ones(size(BF,1)+size(FmBF,1),1),[3]);
  %   Facets.boundary_marker = [ones(size(BF,1),1);-ones(size(FmBF,1),1)];
  %   Facets.holes = cell(numel(Facets.facets),1);
  %   writePOLY_tetgen('temp.poly',V,Facets,[]);
  %   % -p: we're giving PLC, -g: output .mesh, -Y: no steiners
  %   !/usr/local/bin/tetgen -pgY ~/Documents/volume/temp.poly
  %   [VV,TT,FF] = readMESH('temp.1.mesh');
  %
  % Copyright 2011, Alec Jacobson (jacobson@inf.ethz.ch)
  %
  % See also: cdt, tetgen, writePOLY_triangle, writePOLY_pyramid
  %

  v = 1;
  BM = [];
  BMN = [];
  RL = [];
  while v <= numel(varargin)
    switch varargin{v}
    case 'BoundaryMarkers'
      assert((v+1)<=numel(varargin));
      v = v+1;
      BM = varargin{v};
    case 'BoundaryMarkersNodes'
      assert((v+1)<=numel(varargin));
      v = v+1;
      BMN = varargin{v};
    case 'RegionList'
      assert((v+1)<=numel(varargin));
      v = v+1;
      RL = varargin{v};
    otherwise
      error(['Unsupported parameter: ' varargin{v}]);
    end
    v=v+1;
  end

  % open file for writing
  poly_file_handle = fopen(filename,'w');

  % dimensions in V, should be 3
  dim = size(V,2);
  if isempty(V)
    dim = 3;
  end

  if dim ~= 3
    error('writePOLY_tetgen is for 3d meshes. Try writePOLY_triangle etc.');
  end

  % vertices section
  fprintf(poly_file_handle,'# vertices\n');
  fprintf(poly_file_handle,'# Part 1 - node list\n');
  fprintf(poly_file_handle,'%d %d 0 %d\n', size(V,1),size(V,2),~isempty(BMN));
  format = '%d %.17g %.17g %.17g\n';
  if ~isempty(V)
      if ~isempty(BMN)
          assert(numel(BMN)==size(V,1));
          formatV = '%d %.17g %.17g %.17g %d\n';
          fprintf(poly_file_handle,formatV,[1:size(V,1);V';BMN']);
      else
          formatV=format;
          fprintf(poly_file_handle,formatV,[1:size(V,1);V']);

      end
  end

  fprintf(poly_file_handle,'# Part 2 - facet list\n');
  % Try to print all at once if facets are all the same size
  if ~isstruct(F)
    fprintf(poly_file_handle,'%d %d\n',size(F,1),~isempty(BM));
    assert(numel(BM)==size(F,1) || isempty(BM));
    % build format
    fformat = '1 0';
    if ~isempty(BM)
      fformat = [fformat ' %d'];
    end
    fformat = [fformat '\n' num2str(size(F,2))];
    for p=1:size(F,2)
      fformat = [fformat ' %d']; %#ok<AGROW>
    end
    fformat = [fformat '\n'];
    % print all at once
    fprintf(poly_file_handle,fformat,[BM F]');
  else
    fprintf(poly_file_handle,'%d %d\n',numel(F.facets),1);
    % irregular face valences
    for f=1:numel(F.facets)
      % [num polygons] [num holes] [boundary marker]
      fprintf(poly_file_handle,'%d %d %d\n', ...
        numel(F.facets{f}),size(F.holes{f},1),F.boundary_marker(f));
      % loop over polygons
      for p=1:numel(F.facets{f})
        % [num corners] [corner 1] [corner 2] ...
        fprintf(poly_file_handle,' %d',numel(F.facets{f}{p}));
        fprintf(poly_file_handle,' %d',F.facets{f}{p});
        fprintf(poly_file_handle,'\n');
      end
      % [hole #] [hole x] [hole y] [hole z]
      if ~isempty(F.holes{f})
        assert(size(F.holes{f},2) == size(V,2));
        fprintf(poly_file_handle,format,[1:size(F.holes{f},1);F.holes{f}']);
      end
    end
  end

  % [num holes]
  fprintf(poly_file_handle,'# Part 3 - hole list\n');
  fprintf(poly_file_handle,'%d\n',size(H,1));
  if ~isempty(H)
    assert(isempty(V) || size(H,2) == size(V,2));
    fprintf(poly_file_handle,format,[1:size(H,1);H']);
  end
  % [num regions]
  fprintf(poly_file_handle,'# Part 4 - region list\n');
  fprintf(poly_file_handle,'%d\n',size(RL,1));
  if ~isempty(RL)
    % [region #] [region x] [region y] [region z] [region attribute] and/or [region number]
    formatRL = ['%d' repmat(' %.17g',1,size(RL,2)) '\n'];
    fprintf(poly_file_handle,formatRL,[1:size(RL,1);RL']);
  end
  fprintf(poly_file_handle,'\n');
  fclose(poly_file_handle);
end

author: alecjacobson

repo: gptoolbox


text:
%% Winter conference - Fig 3
clear; close all; clc;

Time = {'BN', '24'};
%% MT_textfile load
Path_M = 'H:\1. Journal\NeuroImage\SLEEP_DATA_24H\';
DirGroup_M = dir(fullfile(Path_M,'*'));
FileNamesGroup_M = {DirGroup_M.name};
FileNamesGroup_M = FileNamesGroup_M(1,3:end);
Time_M = {'BN','24'};
removal_word = [".jpg"];
%% load
% Data_BN=importdata('ks_min_20200330_BN_recall_visuo.txt');
ccnt = 0;
for n = 1:size(FileNamesGroup_M,2)
%% PM Task (Succeful)
    % Recall
    for t = 1:size(Time,2)
        Data_recall=textscan(fopen([FileNamesGroup_M{n} '_' Time_M{t} '_recall_visuo.txt']), '%s %s %s %s %s %s %s');
        Data_recall=[Data_recall{:}];
        Hash_Table_recall = str2double(erase(Data_recall([2:end],3),removal_word));
        Hash_Table_recall = [cellfun(@(x) sprintf('%02d',x),num2cell(Hash_Table_recall),'UniformOutput',false)...
            Data_recall([2:end],4) Data_recall([2:end],6) Data_recall([2:end],2)];
        Hash_Table_recall = sortrows(Hash_Table_recall,1);

        for i = 1:size(Hash_Table_recall,1)
            if string(Hash_Table_recall(i,2)) == 'o'
                res(i,t) = 1;
            elseif string(Hash_Table_recall(i,2)) == 'n'
                res(i,t) = 0;
            elseif isspace(string(Hash_Table_recall(i,2))) == 1 % i don't know
                res(i,t) = 2;
            end
        end

        % four possible reponse categories
        temp=[1,0];% 1: old, 0: new
        for i=1:2
            for j=1:2
                if i == 1
                    trial{i,j,t} = find(res(1:38,t)==temp(j));
                elseif i == 2
                    trial{i,j,t} = find(res(39:end,t)==temp(j))+38;
                end
            end
        end
        if t == 1
            Suc = [trial{1,1,t}];
            Succ = sortrows(str2double(Hash_Table_recall(Suc,4)));
        end
    end
%% Picture memory hits
    label = intersect(trial{1,1,1}, trial{1,1,2});

    % Before nap
    Data_recall=textscan(fopen([FileNamesGroup_M{n} '_' Time_M{1} '_recall_visuo.txt']), '%s %s %s %s %s %s %s');
    Data_recall=[Data_recall{:}];
    Hash_Table_recall = str2double(erase(Data_recall([2:end],3),removal_word));
    Hash_Table_recall = [cellfun(@(x) sprintf('%02d',x),num2cell(Hash_Table_recall),'UniformOutput',false)...
        Data_recall([2:end],4) Data_recall([2:end],6) Data_recall([2:end],2)];
    Hash_Table_recall = sortrows(Hash_Table_recall,1);

    Label = sortrows(str2double(Hash_Table_recall(label,4)));
    y = double(ismember(Succ, Label));
    y_ = nonzeros(y);
    STM_hits(n,1) = size(Succ,1);
    LTM_hits(n,1) = size(y_,1);
end

hits = [STM_hits LTM_hits];


text:
function varargout = subsref(S, index)
%SUBSREF    Evaluate a SPINOPERATOR or reference its fields.
%   S(U) is equaivalent to FEVAL(S, U).
%
%   S.PROP returns the property PROP of S.
%
% See also SPINOPERATOR/FEVAL.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

idx = index(1).subs;
switch index(1).type

    case '.'

        % Allow access to any of the properties of S:
        out = S.(idx);

        % Recurse on SUBSREF:
        if ( numel(index) > 1)
            out = subsref(out, index(2:end));
        end

        varargout = {out};

    case '()'

        % Evaluate the operator, i.e., wrapper for SPINOPERATOR/FEVAL:
        varargout = {feval(S, idx{:})};

    otherwise

        error('SPINOPERATOR:subsref:indexType', ...
            ['Unexpected index type encountered: ' index(1).type]);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function f = abs(f)
%ABS   Absolute value of a CHEBFUN3.
%   ABS(F) returns the absolute value of a CHEBFUN3 object F. This function
%   gives an error if the function passes through or becomes numerically
%   close to zero.
%
% See also CHEBFUN3/COMPOSE.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

if ( isempty(f) ) % check for empty CHEBFUN3.
    return
end

if ( isreal(f) )
    % Positive/negative test.
    [ss, ~, ispos] = singleSignTest(f); % Returns TRUE if there is no sign change.
    if ( ss )
        if (ispos)
            return
        else
            f = uminus(f);
            return
        end
    elseif( ~ss )
        error('CHEBFUN:CHEBFUN3:abs:notSmooth', ...
            'Sign change detected. Unable to represent the result.');
    end
else
    % Call the constructor for complex f:
    f = compose(f, @abs);
end

end