text:
function y=unique_and_sort(x,sortopts)
% Unique() function that ignores NaNs in arrays and empty values as well as 'NA' in
% cellstrs, and sorts according to sortopts

%For unique() options
persistent old_matlab_2012
if isempty(old_matlab_2012)
    old_matlab_2012=verLessThan('matlab','7.14');
end

%For iscategorical()
persistent old_matlab_2013
if isempty(old_matlab_2013)
    old_matlab_2013=verLessThan('matlab','8.2');
end

%Create unique in original order
if old_matlab_2012
    [y,ia,~]=unique(x);
    [~,ind]=sort(ia); %Trick to get the original order in older matlab versions
    y=y(ind);
else
    y = unique(x,'stable'); %we keep the original order
end

%Clean up uniques
if ~iscell(x)
    if ~old_matlab_2013 && ~iscategorical(y)
        y(isnan(y)) = []; % remove all nans
    end
else
    y(strcmp(y,'NA'))=[]; %remove all 'NA'
    y(strcmp(y,''))=[]; %remove all ''
end

%Apply sorting options
if numel(sortopts)>1 %Custom ordering
    y=sort(y);%Sort first

    sortopts=shiftdim(sortopts);
    %If correct lengt and we have integers and all numbers from 1 to N are there we probably have indices
    if length(sortopts)==length(y) && isnumeric(sortopts) && sum(sort(sortopts)==(1:length(y))')==numel(y)
        disp('ordering given as indices')
        y=y(sortopts);
    else
        disp('ordering given as values')
        try
            %This should work whatever the lengths of either array
            [present,order]=ismember(sortopts,y);
            y=y(order(present));
        catch
            warning('Improper ordering values')
        end
    end

else %Other orderings
    switch sortopts
        case 1
            y=sort(y); %default
        case -1
            y=flipud(sort(y)); %We use flipud instead of the 'descend' option because somehow it isn't supported for cellstr.
    end
    %case 0 does nothing (keep original ordering)
end

end

author: piermorel

repo: gramm

save_path: github-repos/MATLAB/piermorel-gramm


text: mex gradientMex.cpp

author: flyers

repo: drone-tracking


text:
function HTriggeredRaster(H,V,L);
nPks = 5; % number of examples to plot
pkHt = .1; % min peak height compared to maximum
H(sum(H,2)==0,:)=[];
[N T] = size(V);
[K,~] = size(H);

clf
for fi = 1:K
    axStart = .1 + (fi-1)*.8/K;
    axW = .8/K*.85;
    axL(fi) = subplot('Position', [axStart .1 axW .8]);

    [pks, locs] = findpeaks(H(fi,:),'MinPeakHeight',max(H(fi,:))*pkHt);

    [~, ind] = sort(pks, 'descend');
    ind = ind(1:min(nPks,length(ind)));
    pks = pks(ind); locs = locs(ind);

    R = [];
    twin = [1:L];
    c = 1;
    for pi = 1:length(pks)
        R(c,:,:) = V(:,min(locs(pi)+twin,size(V,2)));
        c = c+1;
    end

    hold on
    for ni = 1:N
        plot([0 L+1], ni*length(pks)*[1 1], 'color', [1 .8 .8], 'linewidth', .1)
    end

    Rmat = reshape(R,N*length(pks), length(twin));
    [Xmat,Ymat] = meshgrid(1:size(Rmat,2),1:size(Rmat,1));
    indBurst = find(Rmat>0);
    s = scatter(Xmat(indBurst)+.9*(-.5+rand(size(indBurst))),Ymat(indBurst), ...
       15, 0*[1 1 1], '.', 'markerfacecolor', 'flat', 'MarkerfaceAlpha',5/nPks, 'MarkerfaceAlpha',5/nPks);

    set(gca, 'ydir', 'reverse')
    axis tight; axis off
end

author: FeeLab

repo: seqNMF

save_path: github-repos/MATLAB/FeeLab-seqNMF


text:
function res = bf_inverse_lcmv(BF, S)
% Computes LCMV filters
% Copyright (C) 2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: bf_inverse_lcmv.m 7703 2019-11-22 12:06:29Z guillaume $

%--------------------------------------------------------------------------
if nargin == 0
    orient        = cfg_menu;
    orient.tag    = 'orient';
    orient.name   = 'Orient to maximum power';
    orient.labels = {'yes', 'no'};
    orient.values = {true, false};
    orient.val    = {true};

    keeplf        = cfg_menu;
    keeplf.tag    = 'keeplf';
    keeplf.name   = 'Keep oriented leadfields';
    keeplf.labels = {'yes', 'no'};
    keeplf.values = {true, false};
    keeplf.val    = {false};

    lcmv      = cfg_branch;
    lcmv.tag  = 'lcmv';
    lcmv.name = 'LCMV';
    lcmv.val  = {orient, keeplf};
    res = lcmv;

    return
elseif nargin < 2
    error('Two input arguments are required');
end

invCy =  BF.features.(S.modality).Cinv;
U     = BF.features.(S.modality).U;

reduce_rank = BF.sources.reduce_rank.(S.modality(1:3));

L = S.L;
W = cell(size(L));

nvert = numel(W);

spm('Pointer', 'Watch');drawnow;
spm_progress_bar('Init', nvert, ['Computing ' S.modality ' filters']); drawnow;
if nvert > 100, Ibar = floor(linspace(1, nvert,100));
else Ibar = 1:nvert; end

for i = 1:nvert
    if ~isnan(L{i})
        lf    = U'*L{i};

        if S.orient
            % Robert's code
            [u, dum] = svd(real(pinv_plus(lf' * invCy *lf, reduce_rank, 0)),'econ');
            eta = u(:,1);
            lf  = lf * eta;
        end

        % construct the spatial filter
        W{i} = pinv_plus(lf' * invCy * lf, reduce_rank, 0)*lf'*invCy;

        if S.keeplf
            L{i} = lf;
        end
    else
        W{i} = NaN;
    end

    if ismember(i, Ibar)
        spm_progress_bar('Set', i); drawnow;
    end
end

spm_progress_bar('Clear');

res.W = W;
res.L = L;

author: spm


text:
function [uOutput] = import_gismo_catalog(nfunction, catalog)
    % import from a GISMO Catalog object
    % Catalog objects are part of the GISMO suite
    %
    % to get helptext:
    %   [helpstring] = import_gismo_catalog(0)
    %
    % to retrieve Catalog stored in file
    %   [uOutput] = import_gismo_catalog(1, filename)
    %
    % to convert Catalog stored in memory
    %   [uOutput] = import_gismo_catalog(1, catalog) % name will be name of variable

    % created by Celso G Reyes, 2017

    % ZMAP format is 10 columns: longitude, latitude, decimal year, month, day,
    % magnitude, depth, hour, minute, second
    uOutput = [];

    if nfunction==0   % Return info about filter
        uOutput = 'GISMO Catalog - import catalog from the GISMO suite';
        return
    end
    if nfunction ==2
        uOutput = 'gismo.html';
        return
    end

    if ~exist('Catalog','class')
        uOutput(['Cannot import catalog, since Catalog class (and probably GISMO) is not installed.',...
            'Last known whereabouts: https://github.com/geoscience-community-codes/GISMO']);
        return
    end
    if ischar(catalog)
        % filename
        if exist(catalog,'file')
            fn = catalog;
            catalog=load(fn);
        end
    end

    if exist('catalog','var') && isa(catalog,'Catalog')
        tb=table(datetime(datevec(catalog.otime)),...
            catalog.lat,...
            catalog.lon,...
            catalog.depth,...
            catalog.mag,...
            catalog.magtype,...
            'VariableNames',{'Date','Latitude','Longitude','Depth','Magnitude','MagnitudeType'});
        uOutput=ZmapCatalog.from(tb);

        uOutput.Name=inputname(2);

    else
        uOutput = 'unable to import Catalog';
    end
end

author: CelsoReyes


text:
% pop_erpparams() - Set plotting and statistics parameters for cluster ERP
%                   plotting
% Usage:
%   >> STUDY = pop_erpparams(STUDY, 'key', 'val');
%
% Inputs:
%   STUDY        - EEGLAB STUDY set
%
% Statistics options:
%   'groupstats  - ['on'|'off'] Compute (or not) statistics across groups.
%                  {default: 'off'}
%   'condstats'  - ['on'|'off'] Compute (or not) statistics across groups.
%                  {default: 'off'}
%   'topotime'   - [real] Plot ERP scalp maps at one specific latency (ms).
%                   A latency range [min max] may also be defined (the
%                   ERP is then averaged over the interval) {default: []}
%   'statistics' - ['param'|'perm'|'bootstrap'] Type of statistics to compute
%                  'param' for parametric (t-test/anova); 'perm' for
%                  permutation-based and 'bootstrap' for bootstrap
%                  {default: 'param'}
%   'naccu'      - [integer] Number of surrogate averages to accumulate for
%                  surrogate statistics. For example, to test whether
%                  p<0.01, use >=200. For p<0.001, use 'naccu' >=2000.
%                  If a threshold (not NaN) is set below, and 'naccu' is
%                  too low, it will be automatically reset. (This option
%                  is now available only from the command line).
%   'threshold'  - [NaN|float<<1] Significance probability threshold.
%                  NaN -> plot the p-values themselves on a different axis.
%                  When possible, the significant time regions are indicated
%                  below the data.
%   'mcorrect'   - ['fdr'|'none'] correction for multiple comparisons
%                  (threshold case only). 'fdr' uses false discovery rate.
%                  See the fdr function for more information. Defaut is
%                  'none'.
% Plot options:
%   'timerange'  - [min max] ERP plotting latency range in ms.
%                  {default: the whole epoch}
%   'ylim'       - [min max] ERP limits in microvolts {default: from data}
%   'plotgroups' - ['together'|'apart'] 'together' -> plot subject groups
%                  on the same axis in different colors, else ('apart')
%                  on different axes. {default: 'apart'}
%   'plotconditions' - ['together'|'apart'] 'together' -> plot conditions
%                  on the same axis in different colors, else ('apart')
%                  on different axes. Note: Keywords 'plotgroups' and
%                  'plotconditions' may not both be set to 'together'.
%                  {default: 'together'}
% See also: std_erpplot()
%
% Authors: Arnaud Delorme, CERCO, CNRS, 2006-

% Copyright (C) Arnaud Delorme, 2006
%
% 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 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

function [ STUDY, com ] = pop_erpparams(STUDY, varargin);

STUDY = default_params(STUDY);
STUDY.etc.erpparams = pop_statparams(STUDY.etc.erpparams, 'default');
TMPSTUDY = STUDY;
com = '';
if isempty(varargin)

    enablecond  = fastif(length(STUDY.design(STUDY.currentdesign).variable(1).value)>1, 'on', 'off');
    enablegroup = fastif(length(STUDY.design(STUDY.currentdesign).variable(2).value)>1, 'on', 'off');
    plotconditions    = fastif(strcmpi(STUDY.etc.erpparams.plotconditions, 'together'), 1, 0);
    plotgroups  = fastif(strcmpi(STUDY.etc.erpparams.plotgroups,'together'), 1, 0);
    vis = fastif(isnan(STUDY.etc.erpparams.topotime), 'off', 'on');

    uilist = { ...
        {'style' 'text'       'string' 'Time range in ms [low high]'} ...
        {'style' 'edit'       'string' num2str(STUDY.etc.erpparams.timerange) 'tag' 'timerange' } ...
        {'style' 'text'       'string' 'Plot limits in uV [low high]'} ...
        {'style' 'edit'       'string' num2str(STUDY.etc.erpparams.ylim) 'tag' 'ylim' } ...
        {'style' 'text'       'string' 'Plot scalp map at latency [ms]' 'enable' vis } ...
        {'style' 'edit'       'string' num2str(STUDY.etc.erpparams.topotime) 'tag' 'topotime' 'enable' vis } ...
        {'style' 'text'       'string' 'Display filter in Hz [high]' } ...
        {'style' 'edit'       'string' num2str(STUDY.etc.erpparams.filter) 'tag' 'filter' } ...
        {} {'style' 'checkbox'   'string' 'Plot first variable on the same panel' 'value' plotconditions 'enable' enablecond  'tag' 'plotconditions' } ...
        {} {'style' 'checkbox'   'string' 'Plot second variable on the same panel' 'value' plotgroups 'enable' enablegroup 'tag' 'plotgroups' } };

    cbline = [0.07 1.1];
    otherline = [ 0.7 .5 0.6 .5];
    geometry = { otherline otherline cbline cbline };

    [STUDY.etc.erpparams res options] = pop_statparams(STUDY.etc.erpparams, 'geometry' , geometry, 'uilist', uilist, ...
                                   'helpcom', 'pophelp(''std_erpparams'')', 'enablegroup', enablegroup, ...
                                   'enablecond', enablecond, ...
                                   'title', 'Set ERP plotting parameters -- pop_erpparams()');

    if isempty(res), return; end;

    % decode inputs
    % -------------
    %if res.plotgroups & res.plotconditions, warndlg2('Both conditions and group cannot be plotted on the same panel'); return; end;
    if res.plotgroups, res.plotgroups = 'together'; else res.plotgroups = 'apart'; end;
    if res.plotconditions , res.plotconditions  = 'together'; else res.plotconditions  = 'apart'; end;
    res.topotime  = str2num( res.topotime );
    res.timerange = str2num( res.timerange );
    res.ylim      = str2num( res.ylim );
    res.filter    = str2num( res.filter );

    % build command call
    % ------------------
    if ~strcmpi( char(res.filter), char(STUDY.etc.erpparams.filter)), options = { options{:} 'filter' res.filter }; end;
    if ~strcmpi( res.plotgroups, STUDY.etc.erpparams.plotgroups), options = { options{:} 'plotgroups' res.plotgroups }; end;
    if ~strcmpi( res.plotconditions , STUDY.etc.erpparams.plotconditions ), options = { options{:} 'plotconditions'  res.plotconditions  }; end;
    if ~isequal(res.ylim     , STUDY.etc.erpparams.ylim),      options = { options{:} 'ylim' res.ylim       }; end;
    if ~isequal(res.timerange, STUDY.etc.erpparams.timerange), options = { options{:} 'timerange' res.timerange }; end;
    if (all(isnan(res.topotime)) & all(~isnan(STUDY.etc.erpparams.topotime))) | ...
            (all(~isnan(res.topotime)) & all(isnan(STUDY.etc.erpparams.topotime))) | ...
                (all(~isnan(res.topotime)) & ~isequal(res.topotime, STUDY.etc.erpparams.topotime))
        options = { options{:} 'topotime' res.topotime };
    end;
    if ~isempty(options)
        STUDY = pop_erpparams(STUDY, options{:});
        com = sprintf('STUDY = pop_erpparams(STUDY, %s);', vararg2str( options ));
    end;
else
    if strcmpi(varargin{1}, 'default')
        STUDY = default_params(STUDY);
    else
        for index = 1:2:length(varargin)
            STUDY.etc.erpparams = setfield(STUDY.etc.erpparams, varargin{index}, varargin{index+1});
        end;
    end;
end;

% scan clusters and channels to remove erpdata info if timerange has changed
% ----------------------------------------------------------
if ~isequal(STUDY.etc.erpparams.timerange, TMPSTUDY.etc.erpparams.timerange)
    if isfield(STUDY.cluster, 'erpdata')
        for index = 1:length(STUDY.cluster)
            STUDY.cluster(index).erpdata  = [];
            STUDY.cluster(index).erptimes = [];
        end;
    end;
    if isfield(STUDY.changrp, 'erpdata')
        for index = 1:length(STUDY.changrp)
            STUDY.changrp(index).erpdata  = [];
            STUDY.changrp(index).erptimes = [];
        end;
    end;
end;

function STUDY = default_params(STUDY)
    if ~isfield(STUDY.etc, 'erpparams'), STUDY.etc.erpparams = []; end;
    if ~isfield(STUDY.etc.erpparams, 'topotime'),    STUDY.etc.erpparams.topotime = []; end;
    if ~isfield(STUDY.etc.erpparams, 'filter'),     STUDY.etc.erpparams.filter = []; end;
    if ~isfield(STUDY.etc.erpparams, 'timerange'),  STUDY.etc.erpparams.timerange = []; end;
    if ~isfield(STUDY.etc.erpparams, 'ylim'     ),  STUDY.etc.erpparams.ylim      = []; end;
    if ~isfield(STUDY.etc.erpparams, 'plotgroups') , STUDY.etc.erpparams.plotgroups = 'apart'; end;
    if ~isfield(STUDY.etc.erpparams, 'plotconditions') ,  STUDY.etc.erpparams.plotconditions  = 'apart'; end;

author: PatternRecognition

repo: OpenBMI


text:
% The COBRAToolbox: testSampleCbModelRHMC.m
%
% Purpose:
%     - tests the sampleCbModel function using RHMC using the E. coli Core Model
%

% initialize the test
fileDir = fileparts(which('testSampleCbModelRHMC'));
models = {};
models{1} = getDistributedModel('ecoli_core_model.mat');
models{2} = getDistributedModel('Acidaminococcus_sp_D21.mat');
load('Recon1.0model.mat', 'Recon1')
models{3} = Recon1;
try
   S = load('iDopaNeuro1.mat');
   names = fieldnames(S);
   models{4} = S.(names{1});
end

fprintf('   Testing sampleCbModel with RHMC ... \n' );

for i = 1:length(models)
   model = models{i};
   n = length(model.lb);
   model.vMean = randn(n,1);
   model.vCov = rand(n,1);
   options.nPointsReturned = 1000;

   for j = 1:3
      [modelSampling, samples, volume] = sampleCbModel(model, '', 'RHMC', options);
      modelSampling.samples = samples;
      ESS = min(effective_sample_size(samples));
      pVal = distribution_test(modelSampling);
      okay = (ESS > 100 && pVal > 0.01 && pVal < 0.99);
      if okay, break; end
   end
   assert(okay);
end

fprintf('Done.\n');

author: opencobra

repo: cobratoolbox


text:
function [PowerPhaseRatemap,spikebinIDs] = bz_PowerPhaseRatemap(spikes,filteredLFP,varargin)
%[PowerPhaseRatemap,spikebinIDs] = bz_PowerPhaseRatemap(spikes,filteredLFP,varargin)
%
%INPUTS
%   spikes          structure containing spikes.times from bz_getSpikes
%   filteredLFP     structure containing filteredLFP.timestamps,
%                       filteredLFP.amp, filteredLFP.phase,
%                       filteredLFP.samplingRate. from bz_Filter(lfp).
%
%  (options)
%   'ints'
%   'powernorm'     (default: 'modZ')
%   'ratenorm'      (default: 'none')
%   'numbins'       (default: 20)
%   'metric'        if you would like to calculate average of some other
%                   metric (i.e. not rate) as function of power/phase
%                   should be cell array (similar to spikes) with value of
%                   other metric for each spike
%   'Nspikethresh'  number of spikes threshold to count mean metric
%                   (default: 10)
%
%OUTPUTS
%   PowerPhaseRatemap
%   spikebinIDs
%
%DLevenstein 2018
%% DEV
p = inputParser;
addParameter(p,'ints',[0 Inf],@isnumeric)
addParameter(p,'powernorm','Zlog')
addParameter(p,'metric',[])
addParameter(p,'Nspikethresh',10)

parse(p,varargin{:})
ints = p.Results.ints;
powernorm = p.Results.powernorm;
metric = p.Results.metric;
Nspikethresh = p.Results.Nspikethresh;

%% Find lfp and spikes in the intervals
instatespiketimes = cellfun(@(X) InIntervals(X,ints),...
    spikes.times,'UniformOutput',false);
instateLFPtimes = InIntervals(filteredLFP.timestamps,ints);

%%
numcells = length(spikes.times);
%%

%Normalize Power
switch powernorm
    case 'Zlog'
        filteredLFP.amp = NormToInt(log10(filteredLFP.amp),'Z',ints,filteredLFP.samplingRate);
end

%%

%Get Power/Phase at each spike
spikes.amp = cellfun(@(X) interp1(filteredLFP.timestamps,filteredLFP.amp,X,'nearest'),...
    spikes.times,'uniformoutput',false);
spikes.phase = cellfun(@(X) interp1(filteredLFP.timestamps,filteredLFP.phase,X,'nearest'),...
    spikes.times,'uniformoutput',false);

%%
numbins = 20;

%Power/Phase Bins
phaseedges = linspace(-pi,pi,numbins+1);
PowerPhaseRatemap.phasebins=phaseedges(1:end-1)+0.5.*diff(phaseedges(1:2));
poweredges = linspace(-1.8,1.8,numbins+1);
PowerPhaseRatemap.powerbins=poweredges(1:end-1)+0.5.*diff(poweredges(1:2));
poweredges(1) = -Inf; poweredges(end) = Inf;

%Calculate Power/Phase Spike Histogram
%phaseamphist = cellfun(@(X,Y) hist3([X,Y],{powerbins,phasebins}),spikeamp,spikephase,'UniformOutput',false);

[PowerPhaseRatemap.Nspikes,~,~,spikebinIDs.powerbin,spikebinIDs.phasebin] = ...
    cellfun(@(X,Y,Z) histcounts2(X(Z),Y(Z),poweredges,phaseedges),...
    spikes.amp,spikes.phase,instatespiketimes,...
    'UniformOutput',false);

%Calculate Power/Phase Time Histogram
%phaseamphist_t = hist3([t_amp,t_phase],{powerbins,phasebins});
%phaseamphist_t = phaseamphist_t./sf_LFP;

PowerPhaseRatemap.occupancy = ...
    histcounts2(filteredLFP.amp(instateLFPtimes),...
    filteredLFP.phase(instateLFPtimes),poweredges,phaseedges);
PowerPhaseRatemap.occupancy = ...
    PowerPhaseRatemap.occupancy./filteredLFP.samplingRate;
PowerPhaseRatemap.powerdist = sum(PowerPhaseRatemap.occupancy,2);

%Normalize Rate
%totaltime = sum(diff(ints,1,2));
%meanrate = cellfun(@(X) length(X)./totaltime,spiketimes,'UniformOutput',false);
PowerPhaseRatemap.ratemap = cellfun(@(X) X./PowerPhaseRatemap.occupancy,...
    PowerPhaseRatemap.Nspikes,'UniformOutput',false);
%phaseamprate = cellfun(@(X,Y) X./Y,phaseamprate,meanrate,'UniformOutput',false);

PowerPhaseRatemap.meanrate = nanmean(cat(3,PowerPhaseRatemap.ratemap{:}),3);

if ~isempty(metric)
    for po = 1:length(PowerPhaseRatemap.powerbins)
        for ph = 1:length(PowerPhaseRatemap.phasebins)

            %Find the spikes in the right bin
            inbinspikes = cellfun(@(X,Y,Z) X==po & Y==ph ,...
                spikebinIDs.powerbin,spikebinIDs.phasebin,...
                'UniformOutput',false);

            %Find the metric from spikes in the state
            instatmetric = cellfun(@(X,Y) X(Y),metric,instatespiketimes,...
                'UniformOutput',false);

            %Map for each cell
            for cc = 1:numcells
                PowerPhaseRatemap.metricmap{cc}(po,ph) = ...
                    nanmean(instatmetric{cc}(inbinspikes{cc}));

                 %Threshold.
                if sum(inbinspikes{cc}) < Nspikethresh
                    PowerPhaseRatemap.metricmap{cc}(po,ph) = nan;
                end
            end

        end
    end

	PowerPhaseRatemap.meanmetric = nanmean(cat(3,PowerPhaseRatemap.metricmap{:}),3);

end
%%
%excell = randsample(spikes.numcells,1);
figure
% subplot(2,2,1)
%     imagesc(PowerPhaseRatemap.phasebins,PowerPhaseRatemap.powerbins,...
%         PowerPhaseRatemap.ratemap{excell})
%     hold on
%     imagesc(PowerPhaseRatemap.phasebins+2*pi,PowerPhaseRatemap.powerbins,...
%         PowerPhaseRatemap.ratemap{excell})
%     xlim([-pi 3*pi])
%     axis xy
%     colorbar

subplot(2,2,1)
    imagesc(PowerPhaseRatemap.phasebins,PowerPhaseRatemap.powerbins,...
        PowerPhaseRatemap.meanrate)
    hold on
    imagesc(PowerPhaseRatemap.phasebins+2*pi,PowerPhaseRatemap.powerbins,...
        PowerPhaseRatemap.meanrate)
    plot(linspace(-pi,3*pi,100),cos(linspace(-pi,3*pi,100)),'k')
    xlim([-pi 3*pi])
    axis xy
    colorbar
    xlabel('Phase');
    ylabel('Power (Z(log))')

subplot(2,2,3)
    bar(PowerPhaseRatemap.powerbins,PowerPhaseRatemap.powerdist)
    xlabel('Power (Z(log))');ylabel('Time (s)')
    box off
    axis tight

if ~isempty(metric)
    subplot(2,2,2)
        h = imagesc(PowerPhaseRatemap.phasebins,PowerPhaseRatemap.powerbins,...
            PowerPhaseRatemap.meanmetric);
        hold on
        h2 = imagesc(PowerPhaseRatemap.phasebins+2*pi,PowerPhaseRatemap.powerbins,...
            PowerPhaseRatemap.meanmetric);
        set(h,'AlphaData',~isnan(PowerPhaseRatemap.meanmetric));
        set(h2,'AlphaData',~isnan(PowerPhaseRatemap.meanmetric));
        plot(linspace(-pi,3*pi,100),cos(linspace(-pi,3*pi,100)),'k')
        xlim([-pi 3*pi])
        axis xy
        colorbar
        xlabel('Phase');
        ylabel('Power (Z(log))')
end

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function test_suite = test_clusterize()
% tests for cosmo_clusterize
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite

function test_clusterize_basics
    ds=cosmo_synthetic_dataset('size','normal','ntargets',1,'nchunks',1);
    sample=ds.samples;

    nh_struct=cosmo_cluster_neighborhood(ds,'progress',false);
    nh=cosmo_convert_neighborhood(nh_struct,'matrix');

    x=sample;
    sample=x>2;
    cl1=cosmo_clusterize(sample,nh);
    assertEqual(cl1,{21,25});
    nb=cosmo_convert_neighborhood(nh,'cell');

    sample=round(x/2);
    cl2=cosmo_clusterize(sample,nh);
    assertEqual(cl2,{[1 11 9 14 17 21 22 23 26 30 25]',...
                    [3 6]',28});

    % test exceptions
    aet=@(varargin)assertExceptionThrown(@()...
                        cosmo_clusterize(varargin{:}),'');
    aet('foo',[]);
    aet(sample,[]);
    aet(sample,-1);
    aet('foo',nh);
    aet(zeros([2 2 2 ]),nh);
    aet([sample;sample],nh);
    aet(sample,zeros([1 1 6]));
    aet(sample,ones(1,7));
    aet(sample,[true;false;true;true;true;false]);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function output = callscs(model)
originalModel = model;
% *********************************************
% Bounded variables converted to constraints
% N.B. Only happens when caller is BNB
% *********************************************
if ~isempty(model.ub)
    [model.F_struc,model.K] = addStructureBounds(model.F_struc,model.K,model.ub,model.lb);
end

% Write nonlinear functions using exponential cone operators, if possible
[model,output] = normalizeExponentialCone(model);
if output.problem
    return
end

% Map to SCS format
data.A = -model.F_struc(:,2:end);
data.b = full(model.F_struc(:,1));
data.c = full(model.c);
cones = [];
cones.f = model.K.f;
cones.l = model.K.l;
cones.q = model.K.q;
cones.s = model.K.s;
cones.ep = model.K.e;

params = model.options.scs;
params.verbose = model.options.verbose;

if params.eliminateequalities
    tempmodel.F_struc = [data.b -data.A];
    tempmodel.c = data.c;
    tempmodel.K = cones;
    [tempmodel,xsol,H] =  removeequalitiesinmodel(tempmodel);
    data.b = full(tempmodel.F_struc(:,1));
    data.A = -tempmodel.F_struc(:,2:end);
    data.c = full(tempmodel.c);
    cones = tempmodel.K;
else
    H = 1;
    xsol = 0;
end

% Extract lower diagonal form for new SCS format
if ~isempty(cones.s) && any(cones.s)
    sdpA = data.A(1+cones.l + cones.f+sum(cones.q):end,:);
    sdpb = data.b(1+cones.l + cones.f+sum(cones.q):end,:);
    expA = data.A(end-3*cones.ep+1:end,:);
    expb = data.b(end-3*cones.ep+1:end,:);
    data.A = data.A(1:cones.l + cones.f+sum(cones.q),:);
    data.b = data.b(1:cones.l + cones.f+sum(cones.q),:);
    top = 1;
    for i = 1:length(cones.s)
        A = sdpA(top:top + cones.s(i)^2-1,:);
        b = sdpb(top:top + cones.s(i)^2-1,:);
        n = cones.s(i);
        ind = find(speye(n));
        b(ind) = b(ind)/sqrt(2);
        A(ind,:) = A(ind,:)/sqrt(2);
        ind = find(tril(ones(n)));
        A = A(ind,:);
        b = b(ind);
        data.A = [data.A;A];
        data.b = [data.b;b];
        top = top  + cones.s(i)^2;
    end
    data.A = [data.A;expA];
    data.b = [data.b;expb];
end

if model.options.savedebug
    save scsdebug data cones params
end

if model.options.showprogress;showprogress(['Calling ' model.solver.tag],model.options.showprogress);end
t = tic;
problem = 0;
switch  model.solver.tag
    case 'scs-direct'
         [x_s,y_s,s,info] = scs_direct(data,cones,params);
    otherwise
        if params.gpu == true
            [x_s,y_s,s,info] = scs_gpu(data,cones,params);
        else
            [x_s,y_s,s,info] = scs_indirect(data,cones,params);
        end
end
solvertime = toc(t);

% Internal format. Only recover the original variables
Primal = H*x_s+xsol;
Primal = Primal(1:length(originalModel.c));
if ~isempty(model.evalMap)
    % No support for duals when exponential cones yet
    Dual = [];
else
    % Map to full format from tril
    Dual = y_s(1:cones.f+cones.l+sum(cones.q));
    if ~isempty(cones.s) && any(cones.s)
        top = 1 + cones.f + cones.l + sum(cones.q);
        for i = 1:length(cones.s)
            n = cones.s(i);
            sdpdual = y_s(top:top + n*(n+1)/2-1,:);
            Z = zeros(n);
            Z(find(tril(ones(n)))) = sdpdual;
            Z = (Z + Z')/2;
            ind = find(speye(n));
            Z(ind) = Z(ind)/sqrt(2);
            Dual = [Dual;Z(:)];
            top = top  + n*(n+1)/2;
        end
    end
end

if nnz(data.c)==0 && isequal(info.status,'Unbounded/Inaccurate')
    info.status = 'Infeasible';
end

switch info.status
    case 'Solved'
        problem = 0;
    case 'Infeasible'
        problem = 1;
     case 'Unbounded'
        problem = 2;
    otherwise
        status = 9;
end

infostr = yalmiperror(problem,model.solver.tag);

% Save ALL data sent to solver
if model.options.savesolverinput
    solverinput.data = data;
    solverinput.cones = cones;
    solverinput.params = params;
else
    solverinput = [];
end

% Save ALL data from the solution?
if model.options.savesolveroutput
    solveroutput.x = x_s;
    solveroutput.y = y_s;
    solveroutput.s = s;
    solveroutput.info = info;
else
    solveroutput = [];
end

% Standard interface
output = createOutputStructure(Primal,Dual,[],problem,infostr,solverinput,solveroutput,solvertime);

function [model,x0,H] = removeequalitiesinmodel(model)

if model.K.f > 0
    % Extract the inequalities
    A_equ = model.F_struc(1:model.K.f,2:end);
    b_equ = -model.F_struc(1:model.K.f,1);

    if 1
        % Find a basis for the column space of A_equ
        [L,U,P] = lu(A_equ');
        r = colspaces(L');
        AA = L';
        H1 = AA(:,r);
        H2 = AA(:,setdiff(1:size(AA,2),r));
        try
            x0 = P'*linsolve(full(L'),linsolve(full(U'),full(b_equ),struct('LT',1==1)),struct('UT',1==1));
        catch
            x0 = A_equ\b_equ;
        end
        % FIX : use L and U stupid!
        H = P'*[-H1\H2;eye(size(H2,2))];

    else
        H = null(full(A_equ));
        x0 = A_equ\b_equ;
    end

    model.c = H'*model.c;
    model.F_struc(:,1) = model.F_struc(:,1) + model.F_struc(:,2:end)*x0;
    model.F_struc = [model.F_struc(:,1) model.F_struc(:,2:end)*H];
    model.F_struc(1:model.K.f,:)=[];
    model.K.f = 0;
else
    H = speye(length(model.c));
    x0 = zeros(length(model.c),1);
end

function  [indx]=colspaces(A)
indx = [];
for i = 1:size(A,2)
    s = max(find(A(:,i)));
    indx = [indx s];
end
indx = unique(indx);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
classdef PnormDescriptor < NormDescriptor

    properties (Access = protected)
        pnorm
    end

    properties (Access = protected, Abstract)
       pos
    end

    methods (Access = protected)

        function computeDistance(obj)
            x = obj.pos;
            p = obj.pnorm;
            xnorm = zeros(size(x,1),1);
            for idim = 1:size(x,2)
               xnorm = xnorm + abs(x(:,idim)).^p;
            end
            d = xnorm.^(1./p);
            obj.dist = d;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%Geometry domain
L=15; H=1; B=1; A=H*B; I = B*H^3/12;
interval = nrbline ([0 0],[L 0]);
problem_data.geo_name = interval;

%Boundary conditions
%Left end: Sliding restraint + concentrated force
%Right end: Simple support + moment
problem_data.drchlt1_ends = [false true]; %deflections, homogeneous
problem_data.drchlt2_ends = [true false]; %rotation angles, homogeneous
problem_data.nmnn1_ends   = [true false]; %forces
problem_data.nmnn2_ends   = [false true]; %moments

M(1) = 0; M(2) = 10;
P = 3;
problem_data.M = M; %moment values
problem_data.P = P; %force value

f_c = 0.5;
problem_data.f = @(x) f_c * sin(pi*x/L); %distributed loading

%Physical parameters
E   = 210000;
problem_data.EI = @(x) E * I* ones (size (x));

% Discretization parameters
method_data.degree     = 3;     % Degree of the basis functions
method_data.regularity = 2;     % Regularity of the basis functions
method_data.nsub       = 20;     % Number of subdivisions
method_data.nquad      = 5;     % Points for the Gaussian quadrature rule

% Solving
[geometry, msh, space, w] = solve_BE_Beam (problem_data, method_data);

%Postprocessing
%Exact solution
c1 = (L*f_c+P*pi)/(pi*E*I);
c2 = -(L^2*f_c+L*P*pi+M(2)*pi)/(pi*E*I);
c3 = -f_c*L^3/(pi^3*E*I);
c4 = L^2*(2*L^2*pi^2*f_c+2*L*P*pi^3+3*M(2)*pi^3+6*L^2*f_c)/(pi^3*E*I*6);

w_ex = @(x) f_c*sin(pi*x/L)/E/I*(L/pi)^4 + c1*x.^3/6 + c2*x.^2/2 + c3*x + c4;
%Gradient of exact solution
gradwex = @(x) f_c*cos(pi*x/L)/E/I*(L/pi)^3 + c1*x.^2/2 + c2*x + c3;
%Hessian of exact solution
hesswex = @(x) -f_c*sin(pi*x/L)/E/I*(L/pi)^2 + c1*x + c2;

%L2, H1 and H2 norms of error

[errh2, errh1, errl2] = sp_h2_error (space, msh, w, w_ex, gradwex, hesswex);
errors = [errh2, errh1, errl2]
[normh2, normh1, norml2] = sp_h2_error (space, msh, zeros (size (w)), w_ex, gradwex, hesswex);
relative_errors = [errh2/normh2, errh1/normh1, errl2/norml2]

[eu, F] = sp_eval (w, space, geometry, 100);
han = figure;
plot(F,eu,F,w_ex(F),'--')
legend('FEM-solution', 'Exact solution')
xlabel('x')
ylabel('w')

%!test
%! ex_BE_beam_static
%! close (han)
%! assert (errors, [1.65638054386650e-06 1.95594407043227e-07 4.14113201313776e-08], 2e-11);
%! assert (relative_errors, [1.31307379889248e-06 1.55064145756187e-07 3.30097917469054e-08], 2e-11);

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [u,perm,gjc,g] = urotorder(u,K, maxu,permIN) %#ok
% [u,perm,gjc,g] = urotorder(u,K, maxu,permIN)
%
% UROTORDER  Stable reORDERing of triu U-factor by Givens ROTations.
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also sedumi

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% 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 should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

%Indicate to the user Matlab cannot find the SeDuMi binaries
sedumi_binary_error();


text:
function w = WeightTransform_Merched(W,E,Dint,N,M)

% WeightTransform_Merched    Mapping of Subband Adaptive Filters Coefficients
%                            to Full-band Filter
%
% Weight trnasformation by Merched, Diniz, and Petraglia (see Section 4.4.3.2)
%
% by Lee, Gan, and Kuo, 2008
% Subband Adaptive Filtering: Theory and Implementation
% Publisher: John Wiley and Sons, Ltd

W = [W, conj(fliplr(W(:,2:end-1)))];
W = W.';
Wp = real(ifft(W,N));

G = zeros(N,M/N);

G(1,:) = Wp(1,1:M/N);                % Discard the last sample
for m = 2:N
    aux = conv(E(m-1,:),Wp(m,:));    % Filtering by subfilters with polyphase
                                     %   components of the analysis filter
    G(m,:) = aux(Dint+2:Dint+1+M/N);
end

w = reshape(G,M,1);	             % Equivalent full-band filter

author: CharlesThaCat

repo: acoustic-interference-cancellation


text:
function z = plus( x, y )
    %PLUS Addition of two TT/MPS operators.
    %   Z = PLUS(X,Y) adds to TT/MPS operators. The resulting TT/MPS operator
    %   has rank equal to the sum of the individual ranks.
    %
    %   TTeMPS Toolbox.
    %   Michael Steinlechner, 2013-2016
    %   Questions and contact: michael.steinlechner@epfl.ch
    %   BSD 2-clause license, see LICENSE.txt

    % add sanity check...
    rx = x.rank;
    ry = y.rank;

    z = TTeMPS_op( cell(1, x.order) );

    % first core:
    tmp = zeros( 1, x.size_col(1), x.size_row(1), rx(2)+ry(2) );
    tmp( 1, :, :, 1:rx(2) ) = x.U{1};
    tmp( 1, :, :, rx(2)+1:end ) = y.U{1};
    z.U{1} = tmp;

    %z.U{1} = reshape( [unfold( x.U{1}, 'left'), unfold( y.U{1}, 'left')], [1, x.size(1), x.rank(2) + y.rank(2)]);

    % central cores:
    for i = 2:x.order-1
        tmp = zeros( rx(i)+ry(i), x.size_col(i), x.size_row(i), rx(i+1)+ry(i+1) );
        tmp( 1:rx(i), :, :, 1:rx(i+1) ) = x.U{i};
        tmp( rx(i)+1:end, :, :, rx(i+1)+1:end ) = y.U{i};
        z.U{i} = tmp;
    end

    % last core:
    tmp = zeros( rx(end-1)+ry(end-1), x.size_col(end), x.size_row(end), 1 );
    tmp( 1:rx(end-1), :, :, 1 ) = x.U{end};
    tmp( rx(end-1)+1:end, :, :, 1 ) = y.U{end};
    z.U{end} = tmp;

    z = update_properties( z );
end

author: NicolasBoumal

repo: manopt


text:
function combo_test01 ( )

%*****************************************************************************80
%
%% COMBO_TEST01 tests BAL_SEQ_ENUM, *_RANK, *_SUCCESSOR, *_UNRANK.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 January 2011
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'COMBO_TEST01\n' );
  fprintf ( 1, '  Balanced sequences:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  BAL_SEQ_ENUM enumerates,\n' );
  fprintf ( 1, '  BAL_SEQ_RANK ranks,\n' );
  fprintf ( 1, '  BAL_SEQ_SUCCESSOR lists,\n' );
  fprintf ( 1, '  BAL_SEQ_UNRANK unranks.\n' );
%
%  Enumerate.
%
  nseq = bal_seq_enum ( n );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  For N = %d\n', n );
  fprintf ( 1, '  the number of balanced sequences is %d\n', nseq );
  fprintf ( 1, '\n' );
%
%  List.
%
  t = [];
  rank = -1;

  while ( 1 )

    rank_old = rank;

    [ t, rank ] = bal_seq_successor ( n, t, rank );

    if ( rank <= rank_old )
      break
    end

    fprintf ( 1, '    %3d  ', rank );
    for i = 1 : 2 * n
       fprintf ( 1, '%2d', t(i) );
    end
    fprintf ( 1, '\n' );

  end
%
%  Unrank.
%
  rank = floor ( nseq / 2 );

  t = bal_seq_unrank ( rank, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Rank = %d\n', rank );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The element of that rank is:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : 2 * n
    fprintf ( 1, '%2d', t(i) );
  end
  fprintf ( 1, '\n' );
%
%  Rank.
%
  rank = bal_seq_rank ( n, t );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The element to be ranked is:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : 2 * n
    fprintf ( 1, '%2d', t(i) );
  end
  fprintf ( 1, '\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Computed rank: %d\n', rank );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [x,state] = struct_gram(z,task)
%STRUCT_GRAM Gramian matrix.
%   [x,state] = struct_gram(z,[]) computes x as the matrix-matrix product
%   z'*z. The structure state stores information which is reused in
%   computing the right and left Jacobian-vector products.
%
%   struct_gram(z,task) computes the right or left Jacobian-vector product
%   of this transformation, depending on the structure task. Use the
%   structure state and add the field 'r' of the same shape as z or the
%   field 'l' of the same shape as x to obtain the structure task for
%   computing the right and left Jacobian-vector products
%
%      (dF(:)/dz(:).')*task.r(:) and
%      (dF(:)/dz(:).')'*task.l(:) + conj((dF(:)/dconj(z(:)).')'*task.l(:)),
%
%   respectively. Here, F(z) represents this transormation, (:) signifies
%   vectorization and the derivative w.r.t. z (conj(z)) is a partial
%   derivative which treats conj(z) (z) as constant. The output has the
%   same shape as x or z for the right and left Jacobian-vector products,
%   respectively.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)
%
%   References:
%   [1] L. Sorber, M. Van Barel, L. De Lathauwer, "Structured data fusion,"
%       ESAT-SISTA Internal Report 13-177, KU Leuven, 2013.

if nargin < 2, task = []; end
state = [];

if isempty(task) || (isempty(task.l) && isempty(task.r))
    x = z'*z;
elseif ~isempty(task.r)
    if ~isreal(z) || ~isreal(task.r)
        error('struct_gram:nonanalytic',['Nonanalytic objective ' ...
            'functions are currently not supported in sdf_nls, please ' ...
            'use sdf_minf instead.']);
    end
    x = z.'*task.r+task.r.'*z;
elseif ~isempty(task.l)
    x = z*(task.l'+task.l);
end

end

author: andrewssobral

repo: lrslibrary


text:
function val=evalBSplinePoly(t,x,k,idx)
%%EVALBSPLINEPOLY Given a set of knots in non-decreasing order (some can be
%                 repeated), evaluate the idx-th order-k B-spline
%                 polynomial. This is B_{idx,k}(x). The first subscript
%                 refers to the knots used in determining the polynomial.
%                 Unlike the function evalBSplinePolys, this function can
%                 handle values of x that are outside of the region
%                 t(idx+-1k)<=x<=t(idx+k).
%
%INPUTS: t A 1XnumKnots or numKnotsX1 vector containing the knots. The
%          knots must be sorted in ascending order. It is possible for the
%          knots to be repeated.
%        x The scalar or matrix set of points at which the polynomials
%          should be evaluated.
%        k Optionally, the order of the B-spline polynomials. If this
%          parameter is omitted, then the default of fix(length(t)/2)+1 is
%          used. This is the highest order that can be used with the given
%          number of knots.
%      idx This input specifies where in terms of the knots the B-spline
%          polynomials are centered. This is an integer value from 1 to
%          numKnots-k.
%
%OUTPUTS: val The value of the b-spline B_{idx,k}(x) evaluated at all of
%             the points in x. This has the same size as x.
%
%This function implements the divided difference formulation of the splines
%from Equation 1 in Chapter IX of [1]. A clearer explanation of the
%notation is given in Chapter 2.4.4 of [2], where Equation 2.4.4.3 is the
%divided difference formula.
%
%EXAMPLE:
%Here, we plot all of the b-spline basis functions from example IX.1 of
%Chapter IX of [1]. The functions only sum to one over the region from 1 to
%6. Thus, the sum is not one over the region plotted before x=1.
% t=[0;1;1;3;4;6;6;6];
% k=3;
% numPoints=500;
% x=linspace(0,6,numPoints);
%
% b1=evalBSplinePoly(t,x,k,1);
% b2=evalBSplinePoly(t,x,k,2);
% b3=evalBSplinePoly(t,x,k,3);
% b4=evalBSplinePoly(t,x,k,4);
% b5=evalBSplinePoly(t,x,k,5);
% figure(1)
% clf
% hold on
% plot(x,b1,'-k','linewidth',2)
% plot(x,b2,'--r','linewidth',2)
% plot(x,b3,'-.g','linewidth',2)
% plot(x,b4,'-b','linewidth',2)
% plot(x,b5,'--c','linewidth',2)
%
%REFERENCES:
%[1] C. de Boor, A Practical Guide to Splines. New York: Springer-Verlag,
%    1978.
%[2] J. Stoer and R. Burlisch, Introduction to Numerical Analysis, 2nd ed.
%    New York: Springer-Verlag, 1991.
%
%April 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

val=zeros(size(x));
numEls=numel(x);
for curIdx=1:numEls
    xCur=x(curIdx);
    tSel=t(idx:(idx+k));
    f=max(tSel-xCur,0).^(k-1);
    fDeriv=@(ti,nd)prod((k-(1:nd)))*max(ti-xCur,0).^(k-1-nd);

    val(curIdx)=(t(idx+k)-t(idx))*evalDividedDiff(tSel,f,fDeriv);
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%Kinematics_String_Replace

r_P1_O_CHAR_X = strrep(r_P1_O_CHAR_X,IN,OUT);
v_P1_O_CHAR_X = strrep(v_P1_O_CHAR_X,IN,OUT);
a_P1_O_CHAR_X = strrep(a_P1_O_CHAR_X,IN,OUT);

r_P2_O_CHAR_X = strrep(r_P2_O_CHAR_X,IN,OUT);
v_P2_O_CHAR_X = strrep(v_P2_O_CHAR_X,IN,OUT);
a_P2_O_CHAR_X = strrep(a_P2_O_CHAR_X,IN,OUT);

r_P1_O_CHAR_Y = strrep(r_P1_O_CHAR_Y,IN,OUT);
v_P1_O_CHAR_Y = strrep(v_P1_O_CHAR_Y,IN,OUT);
a_P1_O_CHAR_Y = strrep(a_P1_O_CHAR_Y,IN,OUT);

r_P2_O_CHAR_Y = strrep(r_P2_O_CHAR_Y,IN,OUT);
v_P2_O_CHAR_Y = strrep(v_P2_O_CHAR_Y,IN,OUT);
a_P2_O_CHAR_Y = strrep(a_P2_O_CHAR_Y,IN,OUT);

author: MatthewPeterKelly

repo: dscTutorials


text:
function p = sigmoid(x)
% Sigmoid function: s(x) = 1 / (1+exp(-x))

% This file is from pmtk3.googlecode.com

p = 1./(1+exp(-x));
end

author: emtiyaz

repo: vadam

save_path: github-repos/MATLAB/emtiyaz-vadam


text:
% The COBRAToolbox: testAddLopLawConstraints.m
%
% Purpose:
%     - test OptForce (according to the tutorial)
%
% Authors:
%     - Thomas Pfau Oct 2017
%

global CBTDIR

solverPkgs = prepareTest('requireOneSolverOf', {'ibm_cplex','gurobi'});

originalDir = pwd;

pathTutorial = which('testOptForce.m');
pathstr = fileparts(pathTutorial);
cd(pathstr)

%Init paralell pool, if possible
try
    minWorkers = 2;
    myCluster = parcluster(parallel.defaultClusterProfile);

    if myCluster.NumWorkers >= minWorkers
        poolobj = gcp('nocreate');  % if no pool, do not create new one.
        if isempty(poolobj)
            parpool(minWorkers);  % launch minWorkers workers
        end
        parPoolCreated = true;
    end
catch
    parPoolCreated = false;
end

model = getDistributedModel('AntCore.mat');
model.c(strcmp(model.rxns, 'R75')) = 1;
model = changeRxnBounds(model, 'EX_gluc', -100, 'l');
model = changeRxnBounds(model, 'EX_o2', -100, 'l');
model = changeRxnBounds(model, 'EX_so4', -100, 'l');
model = changeRxnBounds(model, 'EX_nh3', -100, 'l');
model = changeRxnBounds(model, 'EX_cit', -100, 'l');
model = changeRxnBounds(model, 'EX_glyc', -100, 'l');

origmodel = model;

%set up the xlwrite command for xls io. We do this before the loop, as
%changeCobraSolver will correct the globals which are reset on 2014b by a
%javaaddpath.
setupxlwrite();

for k = 1:length(solverPkgs)
    solverLPOK = changeCobraSolver(solverPkgs.LP{k}, 'LP');
    solverMILPOK = changeCobraSolver(solverPkgs.MILP{k}, 'MILP');

    if solverLPOK && solverMILPOK
        %get max growth rate
        model = origmodel;
        growthRate = optimizeCbModel(model);
        %get max succinate production
        model = changeObjective(model, 'EX_suc');
        maxSucc = optimizeCbModel(model);
        %Set a high Biomass production constraint
        constrWT = struct('rxnList', {{'R75'}}, 'rxnValues', 14, 'rxnBoundType', 'b');
        constrMT = struct('rxnList', {{'R75', 'EX_suc'}}, 'rxnValues', [0, 155.55], ...
                  'rxnBoundType', 'bb');
        succExPos = ismember(model.rxns,'EX_suc');
        biomassPos = ismember(model.rxns,'R75');
        [minFluxesW, maxFluxesW, minFluxesM, maxFluxesM, ~, ~] = FVAOptForce(model, ...
                                                                     constrWT, constrMT);
        %Assert, that the constraints are met.
        assert(minFluxesW(biomassPos) == 14);
        assert(maxFluxesW(biomassPos) == 14);
        assert(minFluxesM(succExPos) == 155.55);
        assert(maxFluxesM(succExPos) == 155.55);
        assert(minFluxesM(biomassPos) == 0);
        assert(maxFluxesM(biomassPos) == 0);
        runID = 'TestOptForceM';

        %Get Must Sets
        fprintf('Building Must sets...\n')
        constrOpt = struct('rxnList', {{'EX_gluc', 'R75', 'EX_suc'}}, 'values', [-100, 0, 155.5]');
        [mustLSet, pos_mustL] = findMustL(model, minFluxesW, maxFluxesW, 'constrOpt', constrOpt, ...
                                  'runID', runID, 'printLevel', 0);
        [mustUSet, pos_mustU] = findMustU(model, minFluxesW, maxFluxesW, 'constrOpt', constrOpt, ...
                                  'runID', runID, 'printLevel', 0);

        %constrOpt = struct('rxnList', {{'EX_gluc', 'R75', 'EX_suc'}}, 'values', [-100, 0, 155.5]');
        exchangeRxns = model.rxns(cellfun(@isempty, strfind(model.rxns, 'EX_')) == 0);
        excludedRxns = unique([mustUSet; mustLSet; exchangeRxns]);
        [mustUU, pos_mustUU, mustUU_linear, pos_mustUU_linear] = ...
                findMustUU(model, minFluxesW, maxFluxesW, 'constrOpt', constrOpt, ...
                           'excludedRxns', excludedRxns,'runID', runID, ...
                           'printLevel', 0);

        [mustLL, pos_mustLL, mustLL_linear, pos_mustLL_linear] = ...
                findMustLL(model, minFluxesW, maxFluxesW, 'constrOpt', constrOpt, ...
                           'excludedRxns', excludedRxns,'runID', runID, ...
                           'printLevel', 0);
        [mustUL, pos_mustUL, mustUL_linear, pos_mustUL_linear] = ...
                findMustUL(model, minFluxesW, maxFluxesW, 'constrOpt', constrOpt, ...
                          'excludedRxns', excludedRxns,'runID', runID, ...
                          'printLevel', 0);
        mustU = unique(union(mustUSet, mustUU));
        mustL = unique(union(mustLSet, mustLL));
        targetRxn = 'EX_suc';
        biomassRxn = 'R75';
        nSets = 1;
        constrOpt = struct('rxnList', {{'EX_gluc','R75'}}, 'values', [-100, 0]);
        fprintf('Running OptForce with k = 1\n')
        [optForceSets, posOptForceSets, typeRegOptForceSets, flux_optForceSets] = ...
                optForce(model, targetRxn, biomassRxn, mustU, mustL, ...
                         minFluxesW, maxFluxesW, minFluxesM, maxFluxesM, ...
                         'k', 1, 'nSets', nSets, 'constrOpt', constrOpt, ...
                         'runID', runID);
        assert(validateOptForceSol(origmodel,posOptForceSets,typeRegOptForceSets,'EX_suc'));
        %clean up
        rmdir(runID,'s')
        fprintf('Running Optforce with k = 2\n')
        nSets = 20;
        runID = 'TestOptForceM2';
        excludedRxns = struct('rxnList', {{'SUCt'}}, 'typeReg','U');
        [optForceSets, posOptForceSets, typeRegOptForceSets, flux_optForceSets] = ...
            optForce(model, targetRxn, biomassRxn, mustU, mustL, ...
             minFluxesW, maxFluxesW, minFluxesM, maxFluxesM, ...
             'k', 2, 'nSets', nSets, 'constrOpt', constrOpt, ...
             'excludedRxns', excludedRxns, ...
             'runID', runID);
        assert(validateOptForceSol(origmodel,posOptForceSets,typeRegOptForceSets,'EX_suc'));
        %clean up
        rmdir(runID,'s');
    end
end

cd(originalDir);


text:
function [P, lam] = pswf(N, c, dom, output_type)
%PSWF   Prolate spheroidal wave functions.
% P = PSWF(N, C) returns a CHEBFUN P representing the Nth prolate spheroidal
% wave function with bandwidth C on the interval [-1,1], i.e., the Nth
% eigenfunction of the differential eigenvalue problem
%
%    [(1-x^2)*P(x)')' + (LAM - C^2*x^2)*P(x) = 0,
%
% with N = 0,1,2,....  C must be a positive scalar but N may be a vector of
% non-negative integers, in which case the output is an array-valued CHEBFUN
% with LENGTH(N) columns. P is scaled so that P'*P = 2/(2N+1), with the sign
% such that sign(P(0)) = (-1)^(N/2) if N is even and sign(P'(0)) =
% (-1)^((N-1)/2) if N is odd (see [3], eq. (30.4.1)).
%
% [P, LAM] = PSWF(N, C) also returns the eigenvalue(s).
%
% Examples:
%    f = pswf(2,pi); f(0.3)
%    plot(pswf(0:2:6, 100))
%
% [1] H. Xiao, V. Rokhlin and N. Yarvin, Prolate spheroidal wavefunctions,
% quadrature and interpolation, Inverse Problems, 17 (2001), 805-838.
%
% [2] https://reference.wolfram.com/language/ref/SpheroidalPS.html
%
% [3] https://dlmf.nist.gov/30.4
%
% See also PSWFPTS, LEGPOLY.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Developer note: The approach is to compute the (approximate) normalised
% Legendre coefficients of the PSWFs by solving an eigenvalue problem [1].
% The Legendre coefficients are then converted to Chebyshev via LEG2CHEB,
% and a Chebfun constructed.
%
% There is functionality for scaling the domain, but this is currently
% undocumented, since we are not sure at present whether this choice
% is the right one.
%
% Copyright 2020 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Defaults:
if ( nargin < 3 )
    dom = [-1 1];
    output_type = 'chebfun';
end
if ( nargin == 3 )
    if ( isnumeric(dom) )
        output_type = 'chebfun';
    else
        output_type = dom;
        dom = [-1 1];
    end
end

% Parse inputs:
assert( nargin >= 2, 'PSWF requires at least two input arguments.')
assert( all(round(N)==N) && all(N>=0) && isvector(N), ...
    'Input N must be vector of non-negative integers.');
assert( (numel(c)==1) && (c>0) , ...
    'Input C must be a positive scalar.');
assert( numel(dom)==2 && all(isfinite(dom)) , ...
    'Domain must be a finite two-vector.');

% Set discretisation size. Heuristic estimates for initialisation.
M = max(ceil([2*sqrt(c)*N, 2*c, 20]));

% Increase discretisation size until the trailing Legendre coefficients are
% sufficiently small:
ishappy = 0;
tol = eps;
count = 0;

while ( ~ishappy )

    % Construct the matrix (see Xiao et al):
    j = (0:M).';
    Asub = c^2*j.*(j-1)./((2*j-1).*sqrt((2*j-3).*(2*j+1)));
    Adia = j.*(j+1) + c^2*(2*j.*(j+1)-1)./((2*j+3).*(2*j-1));
    Asup = c^2*(j+2).*(j+1)./((2*j+3).*sqrt((2*j+5).*(2*j+1)));
    A = diag(Asub(3:end), -2) + diag(Adia, 0) + diag(Asup(1:end-2), 2);

    % Split in to even and odd parts for efficiency/accuracy.
    Ae = A(1:2:end,1:2:end);
    Ao = A(2:2:end,2:2:end);

    % Compute (sorted) eigenvectors:
    [Ve, De] = eig(Ae);
    [lame, idx] = sort(diag(De), 'ascend');
    Ve = Ve(:,idx);
    [Vo, Do] = eig(Ao);
    [lamo, idx] = sort(diag(Do), 'ascend');
    Vo = Vo(:,idx);

    % Reassemble full V and eigenvalues:
    V = zeros(M+1,M+1);
    V(1:2:end,1:2:end) = Ve;
    V(2:2:end,2:2:end) = Vo;
    lam = zeros(M+1,1);
    lam(1:2:end) = lame;
    lam(2:2:end) = lamo;

    % Check discretisation size was large enough;
    ishappy = sum(abs(V(end-3:end,N+1)))/(2*length(N)) < tol;
    if ( ~ishappy )
        M = 2*M;
    end

    % Failsafe:
    count = count + 1;
    if ( count > 10 )
        break
    end

end

% Extract required columns and unnormalise the Legendre coeffients:
V = bsxfun(@times, V(:,N+1), sqrt((0:M)'+1/2) );
lam = lam(N+1);

% Trim trailing coefficients that are below machine precision:
M = max(abs(V), [], 2);
idx = find(M > eps, 1, 'last');
V = V(1:idx,:);

% Scale as per Wolfram Alpha definition [1]:
V = (1./sqrt(N+0.5)).*V;

%%

% Enforce P_N(0) > 0 for N even and P'_N(0) > 0 for N odd.
m = 0:(idx-1)/2; idx = ~mod(N,2);
if ( any(idx) )
    L0 = (-1).^m./(beta(m,.5).*m); L0(1) = 1;
    V0 = L0*V(1:2:end,idx);
    currentSign = sign(V0);
    desiredSign = 1 - mod(N(idx),4);
    scl = currentSign.*desiredSign;
    V(:,idx) = scl.*V(:,idx);
end
if ( ~all(idx) )
    Lp0 = 2*(-1).^m./beta(m+1,.5); Lp0 = Lp0(1:length(V(2:2:end,1)));
    Vp0 = Lp0*V(2:2:end,~idx);
    currentSign = sign(Vp0);
    desiredSign = 1 - mod(N(~idx)-1,4);
    scl = currentSign.*desiredSign;
    V(:,~idx) = scl.*V(:,~idx);
end

% Quit now if only coefficients are required:
if ( strcmpi(output_type, 'coeffs') )
    P = V;
    return
end

% Convert Legendre coeffs to Chebyshev coeffs:
W = leg2cheb(V);

% Enforce even/oddness (which is lost in leg2cheb):
idx = logical(mod(N,2));
W(1:2:end,idx) = 0;
W(2:2:end,~idx) = 0;

% Create a Chebfun from the Chebyshev  coefficients:
P = chebfun(W, dom, 'coeffs');

% The coefficients are trimmed in V, so simplifying should not be necessary.
% P = simplify(P);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function model = read_svm_detection_weights( filepath )
load(filepath, 'weights', 'bias');
model = {weights; bias};
end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
function w = hann(n)
% function w = hann(n)
% A Hann window of length n. Does not smear tones located exactly in a bin.
w = .5*(1 - cos(2*pi*(0:n-1)/n) );

author: Sable

repo: mcbench-benchmarks


text:
function [filename, pathname] = myUiGetFile(browseDir, varargin)
%
% [filename, pathname] = myUiGetFile(browseDir, varargin)
%
% Wrapper for uigetfile. Will start browsing at 'browseDir' instead of pwd.
% The rest of the input args should be the normal arguments to uigetfile (see 'help uigetfile').
%
% HISTORY:
% 2002.03.05 RFD (bob@white.stanford.edu): wrote it.
% 2009.04.29 RAS: allows MultiSelect.
%
% Example:
%
% [fileName,inputGrayPath]=myUiGetFile(get(handles.inputGrayPath,'String'),'*.?ray','Get gray file');
% [fileName,inputGrayPath]=myUiGetFile(inputGrayPath,'*.?ray','Get gray file');
% [fileName,inputGrayPath]=myUiGetFile(inputGrayPath,{'*.?ray','.class'},'Gray or mesh file');
%

curDir = pwd;
if(~exist('browseDir','var')), browseDir = curDir; end

% We're slightly clever- if the browseDir is not a dir, go up one level.
% This lets us pass in a file, and browse the file's dir.
if(~exist(browseDir,'dir'))
    browseDir = fileparts(browseDir);
    if(~exist(browseDir,'dir')), browseDir = pwd; end
end

cd(browseDir);
if(~exist('varargin','var') || isempty(varargin))
    [filename, pathname] = uigetfile('multiselect', 'on');
else
    [filename, pathname] = uigetfile(varargin{:}, 'multiselect', 'on');
end
cd(curDir);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function flat = rmVol2flat(gray,flat)
% flat = rmVol2flat(gray,flat)
%
% sd 04/2008 wrote it

% Don't do this unless gray is really a gray and flat is really a flat
if ~strcmp(gray.viewType,'Gray') || ~strcmp(flat.viewType,'Flat')
    myErrorDlg('vol2flatCorAnal can only be used to transform from gray to flat.');
end

% Check that both gray & flat are properly initialized
if isempty(gray)
  myErrorDlg('Gray view must be open.  Use "Open Gray Window" from the Window menu.');
end
if isempty(flat)
  myErrorDlg('Flat view must be open.  Use "Open Flat Window" from the Window menu.');
end

% Check that dataType is the same for both views. If not, doesn't make sense to do the xform.
% because for example the two dataTypes may have a different number of scans.
checkTypes(gray,flat);

% get/load selected file
try
  rmFile = viewGet(gray,'rmFile');
catch
  disp(sprintf('[%s]:Please select file',mfilename));
  gray    = rmSelect(gray);
  rmFile  = viewGet(gray,'rmFile');
end;
load(rmFile,'model',  'params');
modelGray = model; %#ok<NODEF>

% Mask image for masking the flat map away from where we have data
mask = flat.ui.mask;

% Put up wait bar
waitHandle = mrvWaitbar(0,'Transforming retModel.  Please wait...');
% Intersect the coords from the gray view and the Flat view.
grayIndices=cell(1,2);
flatIndices=cell(1,2);
for h=1:2
    % Get the data corresponding to flat coordinates: coData,
    % ampData, and phData are each of size nVoxels x nScans where
    % nVoxels = size(flat.grayCoords,2).  First, find the
    % intersection of gray.coords and flat.grayCoords.  Then make
    % the data arrays by culling out the values from the
    % intersecting voxels.
    % Find gray nodes that are both in the inplanes and included
    % in the unfold.
    % gray.coords are the gray coords that lie in the inplanes.
    % flat.grayCoords are the gray coords in the unfold.
    % Note: this code segment is essentially identical to
    % code in getFlatCoords.
    [foo,grayIndicesTmp,flatIndicesTmp] = intersectCols(gray.coords,flat.grayCoords{h});
    grayIndices{h}=grayIndicesTmp;
    flatIndices{h}=flatIndicesTmp;
    % Error check on flatIndices.  Because the above code segment
    % is the same as that used to get the flat coords in
    % getFlatCoords, all of the flatIndices should be in the
    % intersection.  If not, something is busted.
    if length(flatIndicesTmp)~=size(flat.grayCoords{h},2)
        myWarnDlg('Ack!  Your flat maps do not appear to come from this segmentation!');
    end
end

ds = viewGet(flat,'datasize');
for m = 1:numel(model),
    % get all model fields
    fnames = fieldnames(model{m});
    for f = 1:numel(fnames),   % loop over model fields
        mrvWaitbar(((m-1).*length(model)+f) ./ (length(model).*length(fnames)))
        paramGray = rmGet(modelGray{m},fnames{f});
        % only xfm model fields that have data of a certain size
        if numel(paramGray)>1 && isnumeric(paramGray),
            % reshape both left and right
            paramFlat = zeros(ds);
            for h=1:2
                % Corresponding coords on the flat map
                coords = flat.coords{h}(:,flatIndices{h});
                switch lower(fnames{f})
                    case {'b',  'beta'}
                        for ii = 1:size(paramGray,1),
                            paramFlat(:,:,h,ii) = ...
								myGriddata(coords, ...
                                           paramGray(ii,grayIndices{h}).',  ...
                                           mask(:,  :,h));
                        end;
                    otherwise,
                        paramFlat(:,  :,h) = ...
							myGriddata(coords, ...
                                       paramGray(grayIndices{h}).', ...
                                       mask(:, :,h));
                end;
            end
            model{m} = rmSet(model{m},fnames{f},paramFlat); %#ok<AGROW>
        end;
    end;
end;
close(waitHandle)

% Now save output
[p outputname] = fileparts(rmFile);
pathStr        = fullfile(dataDir(flat),outputname);

% Overwrite?
if exist(pathStr,'file')
    saveFlag = questdlg([pathStr, ' already exists. Overwrite?'],  ...
						'Save model file?',  'Yes',  'No',  'No');
else
    saveFlag = 'Yes';
end

% Save
if strcmp(saveFlag, 'Yes')
    save(pathStr, 'model', 'params');
    flat = viewSet(flat, 'rmFile', pathStr);
    fprintf(1, '[%s]:Saved %s.\n', mfilename, pathStr);
else
    fprintf(1, '[%s]:Model not saved.\n', mfilename);
end

% If we transformed it, we probably want to view the data. Load the model.
flat = rmSelect(flat, 1, [pathStr '.mat']);
flat = rmLoadDefault(flat);

return;


text:
report_this_filefun(mfilename('fullpath'));

ax = findobj('Tag','main_map_ax');
[x,y, mouse_points_overlay] = select_polygon(ax);
figure_w_normalized_uicontrolunits(map)

plos2 = plot(x,y,'k-','Linewidth',2);        % plot outline
sum3 = 0.;
pause(0.3)

do = [' s' num2str(k) ' = [ x(1) y(1) ; x(2) y(2) ; x(4) y(4) ; x(3) y(3) ]; '];
eval(do)

if k == 1
    do = [' save peruzones.mat  s' num2str(k) ' k  '];
    eval(do)
else
    do = [' save peruzones.mat  s' num2str(k) ' k  -append  '];
    eval(do)
end

cu = cu+1;
do = [ ' s = s' num2str(k) ' ;' ]; eval(do);
p = s; v = s;
% resort polygons from seisrisk format
if length(s(:,1)) == 4;  p(3,:) = s(4,:); p(4,:) = s(3,:); p = [p ; p(1,:) ];  end
if length(s(:,1)) == 6;  p(3,:) = s(4,:); p(4,:) = s(6,:); p(6,:) = s(3,:); p = [p ; p(1,:) ];end
if length(s(:,1)) == 8;  p(3,:) = s(4,:); p(4,:) = s(6,:); p(5,:) = s(8,:); p(6,:) = s(6,:);  p(7,:) = s(5,:); p(8,:) = s(3,:); p = [p ; p(1,:) ];end

pl = plot(p(:,1),p(:,2),'r');
set(pl,'Linewidth',2)

x = [p(:,1)];
y = [p(:,2)];      %  closes polygon

sum3 = 0.;
XI = a.Longitude;          % this substitution just to make equation below simple
YI = a.Latitude;
    l2 = polygon_filter(x,y, XI, YI, 'inside');
newt2 = a.subset(l2);

if newt2.Count > 10   % nur wenn mindestens 6 EQ in zone

    timeplot

    set(pl,'color','k')
    figure_w_normalized_uicontrolunits(map); hold on;
    plot(newt2.Longitude,newt2.Latitude,'go')
    disp(['This is source zone # ' num2str(k) ]);

    %This adjust the data for the completeness and computes rates per year
    l2 = newt2.Date >= 1963; newt3 = newt2(l2,:);td = 30;

    l = newt3(:,6) >= 4.75 & newt3(:,6) < 5.25;
    r0 = length(newt3(l,6))/td;
    l = newt3(:,6) >= 5.25 & newt3(:,6) < 5.75;
    r1 = length(newt3(l,6))/td;
    l = newt3(:,6) >= 5.75 & newt3(:,6) < 6.25;
    r2 = length(newt3(l,6))/td;

    l2 = newt2.Date >= 1930; newt3 = newt2(l2,:); td = 63;

    l = newt3(:,6) >= 6.25 & newt3(:,6) < 6.75;
    r3 = length(newt3(l,6))/td;
    l = newt3(:,6) >= 6.75 & newt3(:,6) < 7.25;
    r4 = length(newt3(l,6))/td;
    l = newt3(:,6) >= 7.25 & newt3(:,6) < 7.95;
    r5 = length(newt3(l,6))/td;

    if r4 == 0 r4 = r5; end
    if r3 == 0 r3 = r4; end
    if r2 == 0 r2 = r3; end
    if r1 == 0 r1 = r2; end
    if r0 == 0 r0 = r1; end

    r = [r0 r1 r2 r3 r4 r5 ];
    f = min(find(r(2:6) == 0));
    %  if isempty(f) == 0
    %   if f < 6  && f > 1
    %      r(f) = r(f-1)*0.0035
    %  end
    %end
    % write info to file
    % cd /home2/stefan/ZMAP/aspar

    s = ['98   1.     -1          zn9930.00' ]; s = s';
    fprintf(fid2,'%s\n',s);
    s = [num2str(length(v)/2,1) ' 1 1']; s = s';
    fprintf(fid2,'%s\n',s);
    s = [num2str(-v(3,1),5) ' ' num2str(-v(3,2),5) ' ' num2str(-v(4,1),5) ' ' num2str(-v(4,2),5)]; s = s';

    fprintf(fid2,'%s\n',s);
    s = [num2str(-v(1,1),5) ' ' num2str(-v(1,2),5) ' ' num2str(-v(2,1),5) ' ' num2str(-v(2,2),5)]; s = s';

    fprintf(fid2,'%s\n',s);

    anz = [r(1) 5.0 ; r(2) 5.5 ; r(3) 6 ; r(4) 6.5 ; r(5) 7.0 ; r(6) 7.5 ];

    fprintf(fid2,'%7.6f    ',anz(:,1));
    fprintf(fid2,'\n');
    fprintf(fid2,'\n');
    fprintf(fid2,'%3.2f  ',anz(:,2));
    fprintf(fid2,'\n');
end %if at leats 6

ans_ = questdlg('  ',...
    'Define More source zones',...
    'Yes please','No thank you','No' );

switch ans_
    case 'Yes please'
        k = k+1;
        defzones
    case 'No thank you'
        return
end

author: CelsoReyes


text:
function scimat2 = scimat_papillary_muscles(scimat, NPAPS)
% SCIMAT_PAPILLARY_MUSCLES  Extract the papillay muscles from a
% segmentation of the Left Ventricle's cavity.
%
% SCIMAT = scimat_papillary_muscles(SCIMAT)
%
%   This function extracts a segmentation of the papillary muscles from the
%   segmentation of a Left Ventricular cavity. In fact, the papillary
%   muscles are extraced, plus some bits of the wall that don't belong to
%   them, and instead of stopping at the chordae tendineae, the
%   segmentation continues and selects some trabeculations.
%
%   SCIMAT is the SCIMAT struct with the segmentation (see "help scimat"
%   for details).
%
% SCIMAT = scimat_papillary_muscles(SCIMAT, NPAPS)
%
%   NPAPS is a constant with the number of papillary muscles we want to
%   extract. By default, NPAPS=2.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2010-2014 University of Oxford
% Version: 0.2.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
narginchk(1, 2);
nargoutchk(0, 1);

% defaults
if (nargin < 2 || isempty(NPAPS))
    NPAPS = 2;
end

% remove the dummy dimension
scimat = scimat_squeeze(scimat);

% get a tight box around the LV segmentation
box = scimat_box(scimat);

% convert real world coordinates to indices
boxi = scimat_world2index(box', scimat.axis)';

% get data volume size
sz = size(scimat.data);

% create output segmentation volume
scimat2 = scimat;
scimat2.data(:) = 0;

%% Compute convex hull for the segmentation in each slice, and remove the
%% segmentation. This way, we get the papillary muscles together with some
%% rubbish where the convex hull fills in incorrectly a region outside the
%% cavity edge

% element for morphological operator
se1 = strel('disk', 1);

% loop every slice with at least one LV cavity voxel
for I = min(boxi(3,:)):max(boxi(3,:))

    % connected components filter (nlabs is the number of labels w/o counting
    % the background)
    [ im, nlabs ] = bwlabel(scimat.data(:, :, I));

    % get number of voxels in each label
    nvox = zeros(nlabs, 1);
    for L = 1:nlabs
        nvox(L) = sum(im(:) == L);
    end

    % we assume that the largest connected component is the cavity
    [ nvoxmax, idx ] = max(nvox);

    % it only makes sense to compute the convex hull if there are at least
    % three pixels
    if (nvoxmax >= 3)

        % we delete the rest of voxels; even if they are part of the cavity, we
        % don't want them interfering with the convex hull
        im = (im == idx);

        % get linear indices of cavity pixels
        idx = find(im);

        % convert linear index to multiple subscripts
        [ir, ic] = ind2sub(sz(1:2), idx);

        % compute the convex hull
        idxhull = convhull(ir, ic);

        % get all the pixels inside of convex hull
        hull = roipoly(im, ic(idxhull), ir(idxhull));

        % we are going to erode the hull polygon 1 pixel, to help in some
        % cases to avoid that the fringe of the next XOR operation are
        % connected between them
        hull = imerode(hull, se1);

        % put the new segmentation into the data volume
        scimat2.data(:, :, I) = bitxor(im, hull);

    end

end

%% Now we are going to get a guess at which bits are papillary muscle and
%% which bits are rubbish. For this, we start start at the middle slice,
%% because this one is likely to show the papillary muscles well. We assume
%% that the largest connected components belong to the papillary muscles,
%% and then look for the components that are touching the previous ones in
%% the next slice

% extract the middle slice
im = scimat2.data(:, :, round(mean(boxi(3,:))));

% compute connected components
[im, nlabs] = bwlabel(im);

% get number of voxels in each label
nvox = zeros(nlabs, 1);
for L = 1:nlabs
    nvox(L) = sum(im(:) == L);
end

% order connected components by decreasing area size
[nvox, idxlabs] = sort(nvox, 1, 'descend');

% the number of papillary muscles at the middle of the cavity may reduce as
% we move towards the end of the volume, and the muscles fuse together
npaps = NPAPS;

% keep only a number of largest components, each one assumed to belong to a
% papillary muscle
im2 = im*0;
for I = 1:npaps
    im2(im == idxlabs(I)) = I;
end

for I = round(mean(boxi(3,:)))+1:max(boxi(3,:))

    % previous slice
    im1 = im2;

    % initialize mask that we are going to use to discard pixels that are
    % far from the pixels in the previous slice
    im1mask = im1*0;

    % dilate each previous papillary muscle component to create a mask for
    % the region of interest
    for L = 1:npaps

        % isolate current component
        im1comp = im1 == L;

        % number of voxels in the component
        nvox = length(find(im1comp));

        % assuming a quasi circular shape, then the area is a=pi*r^2. Because
        % the area is equal to the number of pixels, we can compute the radius
        % of the component as r=sqrt(nvox/pi)
        rad = sqrt(nvox / pi);

        % we are going to allow the component to search for a new component
        % within an area that is the 25% dilation of the old component; but
        % first we do an erosion of 25% to remove spurious bits of
        % segmentation
        im1comp = imerode(im1comp, strel('disk', round(rad*0.25)));

        % dilate the previous component
        im1comp = imdilate(im1comp, strel('disk', round(rad*0.50)));

        % add dilated component to the mask
        im1mask = im1mask | im1comp;
    end

    % remove all pixels from the new slice that are outside of the mask
    scimat2.data(:, :, I) = scimat2.data(:, :, I) & im1mask;

    % compute connected components in the new slice
    [scimat2.data(:, :, I), nlabs2] = bwlabel(scimat2.data(:, :, I));

    % we need a temporal image to store the new papillary bits
    im2 = scimat2.data(:, :, I) * 0;

    % make sure that we are not in a slice without any segmented pixels
    if (nlabs2)

        % matrix to keep track of the intersection of each papillary bit
        % with every connected component
        sim = zeros(npaps, nlabs2);

        % intersect each papillary bit with each connected component bit
        for L = 1:npaps
            % compute how many pixels of the new component overlap the old
            % component
            for C = 1:nlabs2
                sim(L, C) = length(find((im1 == L) ...
                    & (scimat2.data(:, :, I) == C)));
            end
        end

        % for each papillary, select the connected component with the
        % largest intersection with the previous slice's bit
        for L = 1:npaps
            [foo, Lnew] = max(sim(L, :));
            im2(scimat2.data(:, :, I) == Lnew) = L;
        end

    end

    % overrite the segmentation volume with just the papillary bits
    scimat2.data(:, :, I) = im2;

    % recheck the number of components now in the slice (because the
    % papillary muscles can fuse together)
    [ foo, npaps ] = bwlabel(scimat2.data(:, :, I));

end

% keep only a number of largest components, each one assumed to belong to a
% papillary muscle
im2 = im*0;
for I = 1:NPAPS
    im2(im == idxlabs(I)) = I;
end

% the number of papillary muscles at the middle of the cavity may reduce as
% we move towards the end of the volume, and the muscles fuse together
npaps = NPAPS;

for I = round(mean(boxi(3,:))):-1:min(boxi(3,:))

    % previous slice
    im1 = im2;

    % initialize mask that we are going to use to discard pixels that are
    % far from the pixels in the previous slice
    im1mask = im1*0;

    % dilate each previous papillary muscle component to create a mask for
    % the region of interest
    for L = 1:npaps

        % isolate current component
        im1comp = im1 == L;

        % number of voxels in the component
        nvox = length(find(im1comp));

        % assuming a quasi circular shape, then the area is a=pi*r^2. Because
        % the area is equal to the number of pixels, we can compute the radius
        % of the component as r=sqrt(nvox/pi)
        rad = sqrt(nvox / pi);

        % we are going to allow the component to search for a new component
        % within an area that is the 25% dilation of the old component; but
        % first we do an erosion of 25% to remove spurious bits of
        % segmentation
        im1comp = imerode(im1comp, strel('disk', round(rad*0.25)));

        % dilate the previous component
        im1comp = imdilate(im1comp, strel('disk', round(rad*0.50)));

        % add dilated component to the mask
        im1mask = im1mask | im1comp;
    end

    % remove all pixels from the new slice that are outside of the mask
    scimat2.data(:, :, I) = scimat2.data(:, :, I) & im1mask;

    % compute connected components in the new slice
    [scimat2.data(:, :, I), nlabs2] = bwlabel(scimat2.data(:, :, I));

    % we need a temporal image to store the new papillary bits
    im2 = scimat2.data(:, :, I) * 0;

    % make sure that we are not in a slice without any segmented pixels
    if (nlabs2)

        % matrix to keep track of the intersection of each papillary bit
        % with every connected component
        sim = zeros(npaps, nlabs2);

        % intersect each papillary bit with each connected component bit
        for L = 1:npaps
            % compute how many pixels of the new component overlap the old
            % component
            for C = 1:nlabs2
                sim(L, C) = length(find((im1 == L) ...
                    & (scimat2.data(:, :, I) == C)));
            end
        end

        % for each papillary, select the connected component with the
        % largest intersection with the previous slice's bit
        for L = 1:npaps
            [foo, Lnew] = max(sim(L, :));
            im2(scimat2.data(:, :, I) == Lnew) = L;
        end

    end

    % overrite the segmentation volume with just the papillary bits
    scimat2.data(:, :, I) = im2;

    % recheck the number of components now in the slice (because the
    % papillary muscles can fuse together)
    [ foo, npaps ] = bwlabel(scimat2.data(:, :, I));

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function boundary_segment_num = p08_boundary_segment_num ( DUMMY )

%*****************************************************************************80
%
%% P08_BOUNDARY_SEGMENT_NUM counts the boundary segments in problem 08.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Output, integer BOUNDARY_SEGMENT_NUM, the number of boundary segments.
%
  boundary_segment_num = 2;

  return
end


text:
function [tests,pass,perf]=test_wfdbtest(varargin)

inputs={'verbose'};
verbose=0;
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};'])
    end
end
%Sounds redundant, I know, but wfdbtest is used for remote installation
%debugging, it has a  nice suite of test not in the function examples.
%Should be one the first tests run (or for a light check of the toolbox).
test_string={'good=wfdbtest;close all;if(!good) error(''wfdbtest failed'');end'};
[tests,pass,perf]=test_wrapper(test_string,[],verbose);

author: ikarosilva

repo: wfdb-app-toolbox


text:
%MDL_IRB140 Create model of ABB IRB 140 manipulator
%
% MDL_IRB140 is a script that creates the workspace variable irb140 which
% describes the kinematic characteristics of an ABB IRB 140 manipulator
% using standard DH conventions.
%
% Also define the workspace vectors:
%   qz         zero joint angle configuration
%   qr         vertical 'READY' configuration
%   qd         lower arm horizontal as per data sheet
%
% Reference::
% - "IRB 140 data sheet", ABB Robotics.
% - "Utilizing the Functional Work Space Evaluation Tool for Assessing a
%   System Design and Reconfiguration Alternatives"
%   A. Djuric and R. J. Urbanic
%
% Notes::
% - SI units of metres are used.
% - Unlike most other mdl_xxx scripts this one is actually a function that
%   behaves like a script and writes to the global workspace.
%
% See also mdl_fanuc10l, mdl_m16, mdl_motormanHP6, mdl_S4ABB2p8, mdl_puma560, SerialLink.

% MODEL: ABB, IRB140, 6DOF, standard_DH

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function r = mdl_irb140()

    deg = pi/180;

    % robot length values (metres)
    d1 = 0.352;
    a1 = 0.070;
    a2 = 0.360;
    d4 = 0.380;
    d6 = 0.065;

    % DH parameter table
    %     theta d a alpha
    dh = [0 d1 a1  -pi/2
          0 0  a2  0
          0 0  0   pi/2
          0 d4 0   -pi/2
          0 0  0   pi/2
          0 d6 0   pi/2];

    % and build a serial link manipulator

    robot = SerialLink(dh, 'name', 'IRB 140', ...
        'manufacturer', 'ABB', 'ikine', 'nooffset');

    % place the variables into the global workspace
    if nargin == 1
        r = robot;
    elseif nargin == 0
        assignin('caller', 'irb140', robot);
        assignin('caller', 'qz', [0 0 0 0 0 0]); % zero angles
        assignin('caller', 'qd', [0 -90 180 0 0 -90]*deg); % data sheet pose, horizontal
        assignin('caller', 'qr', [0 -90 90 0 90 -90]*deg); % ready pose, arm up
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
%  CUT-N-PASTE
config = configdata(argv()(1));

disp('[Paths]');
disp(strcat('Data: ', config.paths.data));
if isfield(config.paths, 'img')
  disp(strcat('Camera-Images: ', config.paths.img));
end

disp('[Files]');
if isfield(config.files, 'basename')
  disp(strcat('Basename: ', config.files.basename));
end
if isfield(config.files, 'imnames')
  disp(strcat('Image-Name-Prefix: ', config.files.imnames));
end
if isfield(config.files, 'imgext')
  disp(strcat('Image-Extension: ', config.files.imgext));
end

disp('[Images]');
if isfield(config.imgs, 'LEDsize')
  disp(strcat('LED-Size: ', num2str(config.imgs.LEDsize)));
end
if isfield(config.imgs, 'LEDcolor')
  disp(strcat('LED-Color: \'', config.imgs.LEDcolor, '\''));
end
if isfield(config.imgs, 'LEDthr')
  disp(strcat('LED-Threshold: ', num2str(config.imgs.LEDthr)));
end
disp(strcat('Subpix: ', num2str(config.imgs.subpix)));

disp('[Calibration]');
disp(strcat('Nonlinear-Parameters: ', num2str(config.cal.nonlinpar)));
disp(strcat('Nonlinear-Update: ', num2str(config.cal.NL_UPDATE)));
disp(strcat('Initial-Tolerance: ', num2str(config.cal.INL_TOL)));
disp(strcat('Do-Global-Iterations: ', num2str(config.cal.DO_GLOBAL_ITER)));
disp(strcat('Global-Iteration-Threshold: ', num2str(config.cal.GLOBAL_ITER_THR)));
disp(strcat('Global-Iteration-Max: ', num2str(config.cal.GLOBAL_ITER_MAX)));
disp(strcat('Num-Cameras-Fill: ', num2str(config.cal.NUM_CAMS_FILL)));
disp(strcat('Do-Bundle-Adjustment: ', num2str(config.cal.DO_BA)));
disp(strcat('Undo-Radial: ', num2str(config.cal.UNDO_RADIAL)));
disp(strcat('Min-Points-Value: ', num2str(config.cal.MIN_PTS_VAL)));
disp(strcat('N-Tuples: ', num2str(config.cal.NTUPLES)));
disp(strcat('Square-Pixels: ', num2str(config.cal.SQUARE_PIX)));
disp(strcat('Use-Nth-Frame: ',num2str(config.cal.USE_NTH_FRAME)));

author: strawlab

repo: MultiCamSelfCal


text:
function newval=meshinterp(fromval,elemid,elembary,fromelem)
%
% newval=meshinterp(fromval,elemid,elembary,fromelem)
%
% Interpolate nodal values from the source mesh to the target mesh based on
% a linear interpolation
%
% author: Qianqian Fang (q.fang at neu.edu)
%
% input:
%	 fromval: values defined at the source mesh nodes, the row or column
%	          number must be the same as the source mesh node number, which
%	          is the same as the elemid length
%	 elemid: the IDs of the source mesh element that encloses the nodes of
%            the target mesh nodes; a vector of length of target mesh node
%            count; elemid and elembary can be generated by calling
%
%           [elemid,elembary]=tsearchn(node_src, elem_src, node_target);
%
%           note that the mapping here is inverse to that in meshremap()
%
%	 elembary: the bary-centric coordinates of each target mesh nodes
%	         within the source mesh elements, sum of each row is 1, expect
%	         3 or 4 columns (or can be N-D)
%    fromelem: the element list of the source mesh
%
%
% output:
%	 newval: a 2D array with rows equal to the target mesh nodes (nodeto),
%            and columns equals to the value numbers defined at each source
%            mesh node
% example:
%
%    [n1,f1,e1]=meshabox([0 0 0],[10 20 5],1); % target mesh
%    [n2,f2,e2]=meshabox([0 0 0],[10 20 5],2); % src mesh
%    [id, ww]=tsearchn(n2,e2,n1);              % project target to src mesh
%    value_src=n2(:,[2 1 3]);             % create dummy values at src mesh
%    newval=meshinterp(value_src,id, ww, e2);
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if(size(fromval,1)==1)
    fromval=fromval(:);
end

idx=find(~isnan(elemid));

allval=reshape(fromval(fromelem(elemid(idx),:),:),length(idx),size(elembary,2),size(fromval,2));
tmp=cellfun(@(x) sum(elembary(idx,:).*x,2), num2cell(allval,[1 2]),'UniformOutput',false);
newval=size(length(elemid),size(fromval,2));
newval(idx,:)=squeeze(cat(3,tmp{:}));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function Pb=parbound(P,Ps)

% function Pb=parbound(P,Ps)
% ------------------------------------------------------------------------
% Maps the values in the vector P to Pb according to upper and lower limits
% set in the input structure Ps.
%
% Ps.ub -> Upper bound
% Ps.lb -> Lower bound
% Ps.c -> Centre of the output interval (i.e. initial parameter)
% Ps.f -> Scale factor such that at for the input range
% [-f*(ub-lb),f*(ub-lb)] the width of the output range is Ps.t*(ub-lb)
% Ps.t -> Factor (0-1) of the output range width in relation to the intput
% range width
%
%
% %% EXAMPLE
% Ps.f=2;
% Ps.t=0.9;
% Ps.ub=10;
% Ps.lb=0;
% Ps.c=7;
%
% w=Ps.ub-Ps.lb;
% P=linspace(Ps.c-Ps.f*w,Ps.c+Ps.f*w,500);
% Pb=parbound(P,Ps);
% hf1=cFigure; hold on; grid on;
% xlabel('P','FontSize'); ylabel('Pb','FontSize');
% hf=plot(P,Pb,'b-','LineWidth',5);
% hf=plot(P,(Ps.ub-(0.5*(1-Ps.t)*w))*ones(size(P)),'g-','LineWidth',5);
% hf=plot(P,(Ps.lb+(0.5*(1-Ps.t)*w))*ones(size(P)),'g-','LineWidth',5);
% hf=plot(P,Ps.ub*ones(size(P)),'r-','LineWidth',5);
% hf=plot(P,Ps.lb*ones(size(P)),'r-','LineWidth',5);
% axis equal; axis tight; ylim([Ps.lb Ps.ub]);
%
%
%
% Kevin Mattheus Moerman
% kevinmoerman@hotmail.com
% 06/02/2012
%------------------------------------------------------------------------

w=abs(Ps.ub-Ps.lb); %Interval width
b=Ps.f*w; %Point at which the output interval width equals t*w
a=tan(Ps.t*0.5*pi);
x=(P-Ps.c)./(b./a);
Pb=(w.*(0.5+(atan(x)./pi)))+Ps.lb;

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% SP_H1_ERROR: Evaluate the error in H^1 norm.
%
%   [errh1, errl2, errh1s, errh1_elem, errl2_elem, errh1s_elem] = sp_h1_error (space, msh, u, uex, graduex);
%
% INPUT:
%
%     space:    structure representing the space of discrete functions (see sp_scalar/sp_evaluate_col)
%     msh:      structure containing the domain partition and the quadrature rule (see msh_cartesian/msh_evaluate_col)
%     u:        vector of dof weights
%     uex:      function handle to evaluate the exact solution
%     graduex:  function handle to evaluate the gradient of the exact solution
%
% OUTPUT:
%
%     errh1:  error in H^1 norm
%     errl2:  error in L^2 norm
%     errh1s: error in H^1 seminorm
%     errh1_elem:  error in H^1 norm, for each single element
%     errl2_elem:  error in L^2 norm, for each single element
%     errh1s_elem: error in H^1 seminorm, for each single element
%
% Copyright (C) 2010 Carlo de Falco
% Copyright (C) 2011 Rafael Vazquez
% Copyright (C) 2015 Eduardo M. Garau, Rafael Vazquez
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

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

function [errh1, errl2, errh1s, errh1_elem, errl2_elem, errh1s_elem] = sp_h1_error (sp, msh, u, uex, graduex)

  grad_valu = sp_eval_msh (u, sp, msh, 'gradient');
  grad_valu = reshape (grad_valu, sp.ncomp, msh.rdim, msh.nqn, msh.nel);

  for idir = 1:msh.rdim
    x{idir} = reshape (msh.geo_map(idir,:,:), msh.nqn*msh.nel, 1);
  end
  grad_valex  = reshape (feval (graduex, x{:}), sp.ncomp, msh.rdim, msh.nqn, msh.nel);

  w = msh.quad_weights .* msh.jacdet;

  [errl2, errl2_elem] = sp_l2_error (sp, msh, u, uex);
  errh1s_elem = sum (reshape (sum (sum ((grad_valu - grad_valex).^2, 1), 2), [msh.nqn, msh.nel]) .* w);
  errh1s = sqrt (sum (errh1s_elem));

  errh1  = sqrt (errl2^2 + errh1s^2);

  errh1_elem  = sqrt (errl2_elem.^2 + errh1s_elem);
  errh1s_elem = sqrt (errh1s_elem);

end

author: rafavzqz

repo: geopdes


text:
function fh = tapas_physio_plot_cropped_phys_to_acqwindow(ons_secs, sqpar, verbose)
% plot parts of the time series to be processed into regressors
%
% USAGE
%   fh = tapas_physio_plot_cropped_phys_to_acqwindow(ons_secs, sqpar, y)
%
% INPUT
%   ons_secs    - output of tapas_physio_crop_scanphysevents_to_acq_window
%   sqpar       - output of tapas_physio_crop_scanphysevents_to_acq_window
%
% OUTPUT
%   fh          figure handle of output figure

% Author: Lars Kasper
%
% Copyright (C) 2013, Institute for Biomedical Engineering, ETH/Uni Zurich.
%
% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

if nargin == 3
    % If verbose is passed as argument (from updated tapas_physio_review):
    fh = tapas_physio_get_default_fig_params(verbose);
else
    % Backwards compatibility:
    fh = tapas_physio_get_default_fig_params();
end

set(fh,'Name','Preproc: Cutout actual scans - all events and gradients');

Ndummies    = sqpar.Ndummies;
Nslices     = sqpar.Nslices;
sampling    = 1;

%% 1. Plot uncropped data
t        = ons_secs.raw.t;
cpulse      = ons_secs.raw.cpulse;
r           = ons_secs.raw.r;
c           = ons_secs.raw.c;
spulse      = ons_secs.raw.spulse;
svolpulse   = ons_secs.raw.svolpulse;

hasCardiacData = ~isempty(c);
hasRespData = ~isempty(r);

maxValc = max(abs(c));
maxValr = max(abs(r));
if hasCardiacData && hasRespData
    maxVal = maxValc; % max(maxValc, maxValr);
    colors = [1 0 0; 0 1 0];
elseif hasCardiacData
    maxVal = maxValc;
    colors = [1 0 0];

else
    maxVal = maxValr;
    colors = [0 1 0];
end

ampsv = maxVal*1;
amps = maxVal*0.8; % maxVal / 3;
ampc = maxVal*1.2; % maxVal / 2;

%% Plot raw time series data and recorded events as stems
y = [c, r];
x = y(1:sampling:end, :);
stem(spulse(1:Ndummies*Nslices),amps*ones(Ndummies*Nslices,1),'k--');
hold on;
stem(svolpulse(Ndummies+1:end),ampsv*ones(length(svolpulse)-Ndummies,1),'c', 'LineWidth',2);
stem(spulse((Ndummies*Nslices+1):end), amps*ones(length(spulse)-Ndummies*Nslices,1), 'c--') ;

if hasCardiacData
    stem(cpulse, ampc*ones(length(cpulse),1), 'r--') ;
end

for iLine = 1:size(x,2)
   plot(t(1:sampling:end), x(:,iLine), '--', 'Color', colors(iLine,:));
end

%% 2. Plot cropped data

t           = ons_secs.t;
cpulse      = ons_secs.cpulse;
r           = ons_secs.fr;
c           = ons_secs.c;
spulse      = ons_secs.spulse;
svolpulse   = ons_secs.svolpulse;

hs = zeros(1,0);

%plot physiological time courses and scan events
if hasRespData
    hs(end+1) = plot(t,r,'ko');
else
    hs(end+1) = plot(t,c,'ko');
end

y = [c, r];
x = y (1:sampling:end, :);
hs(end+1) = stem(spulse(1:Ndummies*Nslices),amps*ones(Ndummies*Nslices,1),'--k');
hold on;
hs(end+1) = stem(svolpulse(Ndummies+1:end),ampsv*ones(length(svolpulse)-Ndummies,1),'c', 'LineWidth',2);
hs(end+1) = stem(spulse((Ndummies*Nslices+1):end), amps*ones(length(spulse)-Ndummies*Nslices,1), '--c') ;

if hasCardiacData
    hs(end+1) = stem(cpulse, ampc*ones(length(cpulse),1), 'r') ;
end

for iLine = 1:size(x,2)
    hs(end+1) = plot(t(1:sampling:end), x(:,iLine), '-', 'Color', colors(iLine,:));
end

xlabel('t (s)'); ylabel('Amplitude (a. u.)');
title('Cutout region for physiological regressors');

if hasCardiacData && hasRespData

    legend( hs, {
        'used respiratory signal', ...
        ['dummy scan event marker (N = ' int2str(Ndummies*Nslices) ')'], ...
        ['volume event marker (N = ' int2str(length(svolpulse)-Ndummies) '), without dummies'], ...
        ['scan event marker (N = ' int2str(length(spulse)-Ndummies*Nslices) ')'], ...
        ['cardiac pulse (heartbeat) marker (N = ' int2str(length(cpulse)) ')'], ...
        'cardiac signal (dashed = raw)', 'respiratory signal (dashed = raw)'});
elseif hasCardiacData
    legend( hs, {
        'used cardiac signal', ...
        ['dummy scan event marker (N = ' int2str(Ndummies*Nslices) ')'], ...
        ['volume event marker (N = ' int2str(length(svolpulse)-Ndummies) '), without dummies'], ...
        ['scan event marker (N = ' int2str(length(spulse)-Ndummies*Nslices) ')'], ...
        ['cardiac pulse (heartbeat) marker (N = ' int2str(length(cpulse)) ')'], ...
        'cardiac signal (dashed = raw)'});
else % only respData
    legend( hs, {
        'used respiratory signal', ...
        ['dummy scan event marker (N = ' int2str(Ndummies*Nslices) ')'], ...
        ['volume event marker (N = ' int2str(length(svolpulse)-Ndummies) '), without dummies'], ...
        ['scan event marker (N = ' int2str(length(spulse)-Ndummies*Nslices) ')'], ...
        'respiratory signal (dashed = raw)'});
end

ylim(1.4*maxVal*[-1 1]);

author: translationalneuromodeling

repo: tapas


text:
function prob_test078 ( )

%*****************************************************************************80
%
%% TEST078 tests GAMMA_MEAN, GAMMA_SAMPLE, GAMMA_VARIANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  test_num = 2;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST078\n' );
  fprintf ( 1, '  For the Gamma PDF:\n' );
  fprintf ( 1, '  GAMMA_MEAN computes the mean;\n' );
  fprintf ( 1, '  GAMMA_SAMPLE samples;\n' );
  fprintf ( 1, '  GAMMA_VARIANCE computes the variance.\n' );

  a_test(1:2) = [ 1.0, 2.0 ];
  b_test(1:2) = [ 3.0, 0.5 ];
  c_test(1:2) = [ 2.0, 0.5 ];

  for test_i = 1 : test_num

    a = a_test(test_i);
    b = b_test(test_i);
    c = c_test(test_i);

    check = gamma_check ( a, b, c );

    if ( ~check );
      fprintf ( 1, '\n' );
      fprintf ( 1, 'TEST078 - Fatal error!\n' );
      fprintf ( 1, '  The parameters are not legal.\n' );
      return
    end

    mean = gamma_mean ( a, b, c );
    variance = gamma_variance ( a, b, c );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  TEST NUMBER: %6d\n', test_i );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  PDF parameter A =             %14f\n', a );
    fprintf ( 1, '  PDF parameter B =             %14f\n', b );
    fprintf ( 1, '  PDF parameter C =             %14f\n', c );
    fprintf ( 1, '  PDF mean =                    %14f\n', mean );
    fprintf ( 1, '  PDF variance =                %14f\n', variance );

    for i = 1 : nsample
      [ x(i), seed ] = gamma_sample ( a, b, c, seed );
    end

    mean = r8vec_mean ( nsample, x );
    variance = r8vec_variance ( nsample, x );
    xmax = max ( x(1:nsample) );
    xmin = min ( x(1:nsample) );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Sample size =     %6d\n', nsample );
    fprintf ( 1, '  Sample mean =     %14f\n', mean );
    fprintf ( 1, '  Sample variance = %14f\n', variance );
    fprintf ( 1, '  Sample maximum =  %14f\n', xmax );
    fprintf ( 1, '  Sample minimum =  %14f\n', xmin );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function out=isinf(x)

%%%precision=x(1).precision;
%%%out=zeros(size(x));
%%%
%%%for ii=1:numel(x)
%%% imag=false;
%%% [xrval,xival]=getVals(x,ii);
%%% if hasimag(xival), imag=true; end
%%% if imag
%%%  temp=mpfr_isinf(precision,xrval);
%%%  if temp~=0, out(ii)=1; end
%%%  temp=mpfr_isinf(precision,xival);
%%%  if temp~=0, out(ii)=1; end
%%% else
%%%  out(ii)=mpfr_isinf(precision,xrval);
%%%  if out(ii)~=0, out(ii)=1; end
%%% end
%%%end % for ii=1:max(ex,

out=x==inf;

author: opencobra

repo: cobratoolbox