text:
classdef GradientSurfPerimeterComputer < handle

    properties (Access = private)
        outPutFolder
        iMesh
        figureID
        gExperiment
    end

    properties (Access = private)
        inputFiles
        outputFolder
        levelSetParams
        circleCase
    end

    methods (Access = public)

        function obj = GradientSurfPerimeterComputer(cParams)
            obj.init(cParams)
        end

        function compute(obj)
            for im = 1:numel(obj.inputFiles)
                obj.iMesh = im;
                obj.createGradientVariationExperiment();
                nEpsilon = length(obj.gExperiment.regularizedPerimeter.epsilons);
                for iEpsilon = 1:nEpsilon
                    obj.plotSurf(iEpsilon);
                    obj.addXYlabel();
                    obj.addTitle(iEpsilon);
                    obj.printSurf(iEpsilon);
                end
            end
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.inputFiles   = cParams.inputFiles;
            obj.outputFolder = cParams.outputFolder;
            obj.levelSetParams = cParams.levelSetParams;
            obj.circleCase = cParams.circleCase;
        end

        function createGradientVariationExperiment(obj)
            s.levelSetParams = obj.levelSetParams;
            s.inputFile      = obj.inputFiles{obj.iMesh};
            s.iMesh          = obj.iMesh;
            g = GradientVariationExperiment(s);
            obj.gExperiment = g;
        end

        function addXYlabel(obj)
            xlabel('$x$','interpreter','latex')
            ylabel('$y$','interpreter','latex')
        end

        function plotSurf(obj,iEpsilon)
            coord = obj.gExperiment.backgroundMesh.coord;
            x = coord(:,1);
            y = coord(:,2);
            dPer = obj.gExperiment.regularizedPerimeter.perimetersGradient;
            z = dPer(:,iEpsilon);
            tri = delaunay(x,y);
            f = figure();
            h = trisurf(tri,x,y,z);
            shading interp
            obj.figureID = f;
            ar = get(gca,'DataAspectRatio');
            obj.plotBoundaryCutMesh(tri,x,y,z)
            set(gca,'DataAspectRatio',ar);
            view(-60,25);
        end

        function plotBoundaryCutMesh(obj,tri,x,y,z)
            m = obj.computeBoundaryCutMesh();
            %nodes = unique(m.connec(:));
            Xi = m.coord(:,1);
            Yi = m.coord(:,2);
            Zi = interptri(tri,x,y,z,Xi,Yi);
            s.coord = [Xi Yi Zi];
            %s.connec(:,1) = 1:length(Xi);
            %s.connec(:,2) = 2:length(Xi)+1;
            %s.connec(end,2) = 1;
            %s.coord = m.coord;
            s.connec = m.connec;
            s.kFace = -2;
            mB = Mesh(s);
            hold on
            mB.plot
            % axes(axes_h)
            % axis(axis_h)
            % view(a,b)
        end

        function mBCut = computeBoundaryCutMesh(obj)
            g = obj.gExperiment;
            s.backgroundMesh = g.backgroundMesh;
            s.boundaryMesh   = g.boundaryMesh;
            uMesh = UnfittedMesh(s);
            uMesh.compute(g.levelSet);
            bCut = uMesh.boundaryCutMesh;
            mBCut = bCut.mesh;
        end

        function printSurf(obj,iEpsilon)
            part1 = [obj.outputFolder,obj.circleCase,'GradientMesh',num2str(obj.iMesh)];
            part2 = ['Epsilon',num2str(iEpsilon)];
            outFile = [part1,part2];
            printer = surfPrinter(obj.figureID);
            printer.print(outFile)
        end

        function addTitle(obj,iEpsilon)
            epsilons = obj.gExperiment.regularizedPerimeter.epsilons;
            h = obj.gExperiment.backgroundMesh.computeMeanCellSize;
            epsStr = num2str(epsilons(iEpsilon)/h);
            firstStr = 'dPer^R_\varepsilon \ \textrm{with} \ \,';
            secondStr = '\varepsilon/h \, = \, ';
            thirdStr = '\ \textrm{and} \ h/L \, = \, ';
            L = obj.gExperiment.domainLength;
            [n,d] = rat(h/L);
            if n == 1
                hStr = [num2str(n),'/',num2str(d)];
            else
                hStr = num2str(round(h/L,3));
            end
            titleStr = ['$',firstStr,secondStr,epsStr,thirdStr,hStr,'$'];
            title(titleStr,'interpreter','latex');
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function im_ext = bound_extension(im,By,Bx,type);

% im_ext = bound_extension(im,B,type);
%
% Extend an image for avoiding boundary artifacts,
%
%   By, Bx:    widths of the added stripes.
%   type:   'mirror'        Mirror extension
%           'mirror_nr':    Mirror without repeating the last pixel
%           'circular':     fft2-like
%           'zeros'

% Javier Portilla, Universidad de Granada, Jan 2004

[Ny,Nx,Nc] = size(im);

im_ext = zeros(Ny+2*By,Nx+2*Bx,Nc);
im_ext(By+1:Ny+By,Bx+1:Nx+Bx,:) = im;

if strcmp(type,'mirror'),

    im_ext(1:By,:,:) = im_ext(2*By:-1:By+1,:,:);
    im_ext(:,1:Bx,:) = im_ext(:,2*Bx:-1:Bx+1,:);
    im_ext(Ny+1+By:Ny+2*By,:,:) = im_ext(Ny+By:-1:Ny+1,:,:);
    im_ext(:,Nx+1+Bx:Nx+2*Bx,:) = im_ext(:,Nx+Bx:-1:Nx+1,:);
    im_ext(1:By,1:Bx,:) = im_ext(2*By:-1:By+1,2*Bx:-1:Bx+1,:);
    im_ext(Ny+1+By:Ny+2*By,Nx+1+Bx:Nx+2*Bx,:) = im_ext(Ny+By:-1:Ny+1,Nx+Bx:-1:Nx+1,:);
    im_ext(1:By,Nx+1+Bx:Nx+2*Bx,:) = im_ext(2*By:-1:By+1,Nx+Bx:-1:Nx+1,:);
    im_ext(Ny+1+By:Ny+2*By,1:Bx,:) = im_ext(Ny+By:-1:Ny+1,2*Bx:-1:Bx+1,:);

elseif strcmp(type,'mirror_nr'),

    im_ext(1:By,:,:) = im_ext(2*By+1:-1:By+2,:,:);
    im_ext(:,1:Bx,:) = im_ext(:,2*Bx+1:-1:Bx+2,:);
    im_ext(Ny+1+By:Ny+2*By,:,:) = im_ext(Ny+By-1:-1:Ny,:,:);
    im_ext(:,Nx+1+Bx:Nx+2*Bx,:) = im_ext(:,Nx+Bx-1:-1:Nx,:);
    im_ext(1:By,1:Bx,:) = im_ext(2*By+1:-1:By+2,2*Bx+1:-1:Bx+2,:);
    im_ext(Ny+1+By:Ny+2*By,Nx+1+Bx:Nx+2*Bx,:) = im_ext(Ny+By-1:-1:Ny,Nx+Bx-1:-1:Nx,:);
    im_ext(1:By,Nx+1+Bx:Nx+2*Bx,:) = im_ext(2*By+1:-1:By+2,Nx+Bx-1:-1:Nx,:);
    im_ext(Ny+1+By:Ny+2*By,1:Bx,:) = im_ext(Ny+By-1:-1:Ny,2*Bx+1:-1:Bx+2,:);

elseif strcmp(type,'circular'),

    im_ext(1:By,:,:) =  im_ext(Ny+1:Ny+By,:,:);
    im_ext(:,1:Bx,:) = im_ext(:,Nx+1:Nx+Bx,:);
    im_ext(Ny+1+By:Ny+2*By,:,:) = im_ext(By+1:2*By,:,:);
    im_ext(:,Nx+1+Bx:Nx+2*Bx,:) = im_ext(:,Bx+1:2*Bx,:);
    im_ext(1:By,1:Bx,:) = im_ext(Ny+1:Ny+By,Nx+1:Nx+Bx,:);
    im_ext(Ny+1+By:Ny+2*By,Nx+1+Bx:Nx+2*Bx,:) = im_ext(By+1:2*By,Bx+1:2*Bx,:);
    im_ext(1:By,Nx+1+Bx:Nx+2*Bx,:) = im_ext(Ny+1:Ny+By,Bx+1:2*Bx,:);
    im_ext(Ny+1+By:Ny+2*By,1:Bx,:) = im_ext(By+1:2*By,Nx+1:Nx+Bx,:);

end

author: ricedsp


text:
% Function to derive residual peak prominence parameter
%
% Description
% Function to derive residual peak prominence parameter
%
%
% Inputs
%  x        : [samples] [Nx1]  Speech signal
%  fs       : [Hz]      [1x1]  Sampling frequency
%  res      : [samples] [Nx1]  Linear prediction residual of speech
%  Es       : [dB]      [Mx1]  Energy contour
%
% Outputs
%  peak_inter:[samples] [Nx2]  Interpolated peak prominence contour
%
% Example
%  Please see the HOWTO_glottalsource.m example file.
%
% References
%  [1] Drugman, T., Kane, J., Gobl, C., `Automatic Analysis of Creaky
%       Excitation Patterns', Submitted to Computer Speech and
%       Language.
%  [2] Kane, J., Drugman, T., Gobl, C., (2013) `Improved automatic
%       detection of creak', Computer Speech and Language 27(4), pp.
%       1028-1047.
%  [3] Drugman, T., Kane, J., Gobl, C., (2012) `Resonator-based creaky
%       voice detection', Interspeech 2012, Portland, Oregon, USA.
%
% Copyright (c) 2013 University of Mons, FNRS & 2013 Trinity College Dublin
%
% License
%  This code is a part of the GLOAT toolbox with the following
%  licence:
%  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.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Authors
%  Thomas Drugman <thomas.drugman@umons.ac.be> & John Kane <kanejo@tcd.ie>

function [peak_inter,peak_prom,peak_t,rep] = res_peak(x,fs,F0mean,res,Es)

% Function to generate residual peak prominence contour which makes up the
% second component of the creak detection algorithm.

% Different settings according to the speakers baseline pitch
if F0mean>100
    maxCreakF0=90;
elseif F0mean < 100 && F0mean>=85
    maxCreakF0=65;
elseif F0mean<85
    maxCreakF0=55;
else maxCreakF0=80;
end

% Set window length based on maximum possible creaky F0
winLen=round(fs/maxCreakF0)*2;

% Resonator settings
Phi=2*pi*1*F0mean/fs;
Rho=0.8;
rep=filter([1 0 0],[1 -2*Rho*cos(Phi) Rho^2],res);

% Measure residual peak prominence
[peak_prom,peak_t] = get_res_peak_prom(rep,fs,winLen,x,Es);

% Interpolate
if length(peak_prom)>1
    peak_inter=interp1(peak_t,peak_prom,1:length(x));
else peak_inter=zeros(1,length(x));
end

peak_inter=[peak_inter(:) (1:length(x))'];

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function result = outOfChina(lat,long)

    result = [long < 72.004 | long > 137.8347 | lat < 0.8293 | lat > 55.8271];

end

author: googollee

repo: eviltransform


text:
%% DOS_DIR_TEST tests the DOS facility for issuing operating system commands.
%
%  Discussion:
%
%    DIR is a legal command on MS/DOS systems, and returns a list of the
%    files in the current directory.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 July 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'DOS_DIR_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Demonstrate a simple use of the "dos" command,\n' );
  fprintf ( 1, '  which allows MATLAB to invoke MS/DOS operating system commands.\n' );

  [ status, string ] = dos ( 'dir' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Command status was %d\n', status );
  fprintf ( 1, '  Command output was %s\n', string );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Repeat the command, but add the -ECHO option.\n' );

  [ status, string ] = dos ( 'dir', '-echo' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Command status was %d\n', status );
  fprintf ( 1, '  Command output was %s\n', string );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now repeat command, but with a final ampersand.\n' );

  [ status, string ] = dos ( 'dir &' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Command status was %d\n', status );
  fprintf ( 1, '  Command output was %s\n', string );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'DOS_DIR_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

author: johannesgerer

repo: jburkardt-m


text:
function test_suite=test_meeg_source()
% tests for MEEG datasets in source space
%
% #   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_meeg_dataset()
    tps={'freq_pow','time_mom','rpt_trial_mom','rpt_trial_mom_1d',...
                        'time_mom_1d',...
                        'rpt_trial_pow','strange_ft_lcmv_pow'};
    for j=1:numel(tps)
        is_ft_strange_lcmv_avg_pow=isequal(cosmo_strsplit(tps{j},'_',1),...
                                                'strange');

        [ft,fdim,data_label]=generate_ft_source(tps{j});
        ds=cosmo_meeg_dataset(ft);

        % check fdim
        assertEqual(ds.a.fdim,fdim);

        key=data_label{1};
        sub_key=data_label{2};

        ft_first_sample_all=ft.(key).(sub_key);

        if is_ft_strange_lcmv_avg_pow
            ft_first_sample=ft_first_sample_all(:,ft.inside);
        else
            ft_first_sample=ft_first_sample_all(ft.inside,:);
        end

        switch sub_key
            case 'mom'
                % choose a random sensor
                inside_idxs=find(ft.inside);
                pos=ceil(rand()*numel(inside_idxs));
                inside_idx=inside_idxs(pos);
                assertEqual(ds.samples(1,ds.fa.pos==inside_idx),...
                                        ft_first_sample{pos}(:)')

            otherwise
                assertEqual(ds.samples(1,:),ft_first_sample(:)');
        end

        [ft_arr, ft_labels, ft_values]=cosmo_unflatten(ds,2,...
                                            'matrix_labels','pos');

        if is_ft_strange_lcmv_avg_pow
            ft_arr=ft_arr';
        end
        assertEqual(ft_arr(ft_arr~=0),ds.samples(:));
        assertEqual(ft_labels, fdim.labels);
        assertEqual(ft_values, fdim.values);

        % select single element, and ensure it is the same in the
        % fieldtrip struct as in the dataset struct
        dim_sizes=cellfun(@(x)size(x,2),fdim.values);
        ndim=numel(dim_sizes);
        rp=ceil(rand(1,ndim).*dim_sizes(:)');
        [nsamples,nfeatures]=size(ds.samples);
        ds_msk=false(1,nfeatures);
        ft_idx=cell(1,1+ndim);
        ft_idx{1}=randperm(nsamples);
        for k=1:ndim
            dim_label=fdim.labels{k};
            ds_msk = ds_msk | rp(k)~=ds.fa.(dim_label);

            switch dim_label
                case 'mom'
                    ft_idx{k+1}=rp(k);

                case 'pos'
                    ft_values=ft.(dim_label);
                    ft_idx{k+1}=find(all(bsxfun(@eq,...
                                    ft_values(rp(k),:),ft_values),2));

                otherwise
                    ft_values=ft.(dim_label);
                    ft_idx{k+1}=find(ft_values(rp(k))==ft_values);
            end
        end

        ds_sel=cosmo_slice(cosmo_slice(ds,~ds_msk,2),ft_idx{1});

        if is_ft_strange_lcmv_avg_pow
            ft_idx=ft_idx(2:end);
        end

        ft_sel=ft_arr(ft_idx{:});

        if ft_sel==0
            assertTrue(isempty(ds_sel.samples));
        else
            assertEqual(ds_sel.samples,ft_sel);
        end

        %re-order features
        nfeatures=size(ds.samples,2);
        ds2=cosmo_slice(ds,randperm(nfeatures),2);

        ft2=cosmo_map2meeg(ds2);

        if is_ft_strange_lcmv_avg_pow
            ft2.avg.pow=ft2.avg.pow';
            ft2.time=ft.time;
        end

        assertEqual(ft,ft2);

        if j==1
            % test compatibility with old fieldtrip

            ft2.inside=find(ft2.inside);
            ds3=cosmo_meeg_dataset(ft2);
            assertEqual(ds,ds3);

            ft2.inside=struct();
            assertExceptionThrown(@()cosmo_meeg_dataset(ft2),'');
        end
    end

function test_meeg_fmri_dataset()
    ds=cosmo_synthetic_dataset('type','source');
    ds_fmri=cosmo_fmri_dataset(ds);
    ft=cosmo_map2meeg(ds);
    ds_ft_fmri=cosmo_fmri_dataset(ft);

    ds_vol=cosmo_vol_grid_convert(ds,'tovol');
    assertEqual(ds_vol,ds_fmri);

    assertTrue(isempty(fieldnames(ds_ft_fmri.sa)))
    ds_vol=rmfield(ds_vol,'sa');
    ds_ft_fmri=rmfield(ds_ft_fmri,'sa');
    assertEqual(ds_vol,ds_ft_fmri);

function test_irregular_source_grid()
    ft=generate_ft_source('freq_pow');
    ft.pos=ft.pos+randn(size(ft.pos));

    ds=cosmo_meeg_dataset(ft);
    ft2=cosmo_map2meeg(ds);

    assertEqual(ft,ft2);

function [ft,fdim,data_label]=generate_ft_source(tp)
    ft=struct();
    dim_pos_range={-3:3,-4:4,-5:5};
    nsamples=2;
    freq=[3 5 7 9];
    time=[-1 0 1 2];
    mom_labels_3d={'x','y','z'};
    mom_labels_1d={'xyz'};

    ft.dim=cellfun(@numel,dim_pos_range);
    ft.pos=cosmo_cartprod(dim_pos_range);
    ft.inside=sum(ft.pos.^2,2)<30;

    fdim=struct();

    switch tp
        case 'freq_pow'
            ft.freq=freq;
            ft.method='average';
            ft.avg=generate_data(ft.inside,numel(freq),1,'pow');
            fdim.labels={'pos';'freq'};
            fdim.values={ft.pos';freq(:)'};
            data_label={'avg','pow'};

        case 'time_mom'
            ft.time=time;
            ft.method='average';
            ft.avg=generate_data(ft.inside,numel(time),1,'mom3d');
            fdim.labels={'pos';'mom';'time'};
            fdim.values={ft.pos';mom_labels_3d;time(:)'};
            data_label={'avg','mom'};

        case 'time_mom_1d';
            ft.time=time;
            ft.method='average';
            ft.avg=generate_data(ft.inside,numel(time),1,'mom1d');
            fdim.labels={'pos';'mom';'time'};
            fdim.values={ft.pos';mom_labels_1d;time(:)'};
            data_label={'avg','mom'};

        case 'rpt_trial_pow';
            ft.time=time;
            ft.method='rawtrial';
            ft.trial=generate_data(ft.inside,numel(time),nsamples,'pow');
            fdim.labels={'pos';'time'};
            fdim.values={ft.pos';time(:)'};
            data_label={'trial','pow'};

        case 'rpt_trial_mom';
            ft.time=time;
            ft.method='rawtrial';
            ft.trial=generate_data(ft.inside,numel(time),nsamples,'mom3d');
            fdim.labels={'pos';'mom';'time'};
            fdim.values={ft.pos';mom_labels_3d;time(:)'};
            data_label={'trial','mom'};

        case 'rpt_trial_mom_1d';
            ft.time=time;
            ft.method='rawtrial';
            ft.trial=generate_data(ft.inside,numel(time),nsamples,'mom1d');
            fdim.labels={'pos';'mom';'time'};
            fdim.values={ft.pos';mom_labels_1d;time(:)'};
            data_label={'trial','mom'};

        case 'strange_ft_lcmv_pow'
            % FT with LCMV source puts .avg data in a 1xNSOURCE matrix,
            % even with NTIME>1 timepoints. In this case, the output
            % dataset structure has the time points averaged.
            ft.time=time;
            ft.method='average';
            ft.avg=generate_data(ft.inside,1,1,'pow');
            ft.avg.pow=ft.avg.pow';
            fdim.labels={'pos';'time'};
            fdim.values={ft.pos';mean(time)};
            data_label={'avg','pow'};

        otherwise
            error('unsupported type %s', tp);

    end

function d=generate_data(inside,nfeatures,nsamples,fld)
    is_single_trial=nargin<3;

    if is_single_trial
        nsamples=1;
    end

    switch fld
        case 'mom3d'
            d=generate_mom(inside,nfeatures,nsamples,3);
        case 'mom1d'
            d=generate_mom(inside,nfeatures,nsamples,1);
        case 'pow'
            d=generate_pow(inside,nfeatures,nsamples);
        otherwise
            error('not supported: %s', fld);
    end

function all_trials=generate_pow(inside,nfeatures,nsamples)
    nf=numel(inside);
    ni=sum(inside);

    all_data=NaN(nf,nfeatures);
    trial_data=cosmo_rand(ni,nfeatures);

    all_trials_cell=cell(nsamples,1);
    for j=1:nsamples
        data=all_data;
        data(inside,:)=trial_data+j;
        all_trials_cell{j}.pow=data;
    end

    all_trials=cat(1,all_trials_cell{:})';

function all_trials=generate_mom(inside,nfeatures,nsamples,ndim)
    nf=numel(inside);
    i=find(inside);
    ni=numel(i);

    data=cosmo_rand(ndim,nfeatures,ni);

    all_trials_cell=cell(1,nsamples);
    for j=1:nsamples
        one_trial=cell(nf,1);
        for k=1:ni
            one_trial{i(k)}=data(:,:,k)+j;
        end
        all_trials_cell{j}=one_trial;
    end

    all_trials=struct('mom',all_trials_cell);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [signal, lineNoiseOut] = removeLineNoise(signal, lineNoiseIn)
% Remove sharp spectral peaks from signal using Sleppian filters
%
% Usage:
% signal = cleanLineNoise(signal)
% [signal, lineNoiseOut] = hcleanLineNoise(signal, lineNoiseIn)
%
% Parameters:
%    signal          Structure with .data and .srate fields
%    lineNoiseIn     Input structure with fields described below
%
% Structure parameters (lineNoiseIn):
%    fPassBand       Frequency band used (default [0, Fs/2] = entire band)
%    Fs 	            Sampling frequency
%    fScanBandWidth  +/- bandwidth centered on each f0 to scan for significant
%                       lines (TM)
%    lineFrequencies Line frequencies to be removed (default
%                       [60, 120, 180, 240, 300])
%    lineNoiseChannels  Channels to remove line noise from (default
%                       size(data, 1))
%    maximumIterations   Maximum times to iterate removal (default = 10)
%    p               Significance level cutoff (default = 0.01)
%    pad             FFT padding factor ( -1 corresponds to no padding,
%                       0 corresponds to padding to next highest power of 2
%                       etc.) (default is 0)
%    pnts
%    tapers          Precomputed tapers from dpss
%    taperBandWidth  Taper bandwidth (default 2 Hz)
%    taperWindowSize Taper sliding window length (default 4 sec)
%    taperWindowStep Sliding window step size (default 4 sec = no overlap)
%    tau             Window overlap smoothing factor (default 100)
%
% This function is based on code originally written by Tim Mullen in a
% package called tmullen-cleanline which is based on the chronux_2
% libraries.
%
%% Check the incoming parameters
if nargin < 1
    error('removeLineNoise:NotEnoughArguments', 'requires at least 1 argument');
elseif isstruct(signal) && ~isfield(signal, 'data')
    error('removeLineNoise:NoDataField', 'requires a structure data field');
elseif size(signal.data, 3) ~= 1
    error('removeLineNoise:DataNotContinuous', 'signal data must be a 2D array');
elseif size(signal.data, 2) < 2
    error('removeLineNoise:NoData', 'signal data must have multiple points');
elseif ~exist('lineNoiseIn', 'var') || isempty(lineNoiseIn)
    lineNoiseIn = struct();
elseif isempty(lineNoiseIn) || ~isstruct(lineNoiseIn)
    error('removeLineNoise:NoData', 'second argument must be a structure')
end

%% Set the defaults to appropriate values
defaults = getPrepDefaults(signal, 'linenoise');
lineNoiseOut = struct('lineNoiseMethod', [], ...
    'lineNoiseChannels', [], 'Fs', [], ...
    'lineFrequencies', [], 'p', [], 'fScanBandWidth', [], ...
    'taperBandWidth', [], 'taperWindowSize', [], ...
    'taperWindowStep', [], 'tau', [], 'pad', [], ...
    'fPassBand', [], 'maximumIterations', []);

[lineNoiseOut, errors] = checkPrepDefaults(lineNoiseIn, lineNoiseOut, defaults);
if ~isempty(errors)
    error('removeLineNoise:BadParameters', ['|' sprintf('%s|', errors{:})]);
end

if strcmpi(lineNoiseOut.lineNoiseMethod, 'clean')
    [signal, lineNoiseOut] = cleanLineNoise(signal, lineNoiseOut);
elseif strcmpi(lineNoiseOut.lineNoiseMethod, 'blasst')
    [signal, lineNoiseOut] = blasstLineNoise(signal, lineNoiseOut);
elseif ~strcmpi(lineNoiseOut.lineNoiseMethod, 'none')
    error('removeLineNoise:BadLineNoiseMethod', ...
          'Unrecognized line noise removal method');
end

author: VisLab

repo: EEG-Clean-Tools


text:
%%***************************************************************************
%% steplength: compute xstep such that  X + xstep*dX >= 0.
%%
%% [xstep] = steplength(blk,X,dX,Xchol,invXchol);
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%***************************************************************************

    function [xstep,invXchol] = steplength(blk,X,dX,Xchol,invXchol);

%%
    for p = 1:size(blk,1)
       pblk = blk(p,:);
       numblk = length(pblk{2});
       pblksize = sum(pblk{2});
       if (any(isnan(dX{p})) | any(isinf(dX{p}))); xstep = 0; break; end;
       if strcmp(pblk{1},'s')
          if (max(pblk{2}) >= 200)
             use_lanczos = 1;
          else
             use_lanczos = 0;
          end
          if (use_lanczos)
             tol = 1e-3;
             maxit = max(min(pblksize,30),round(sqrt(pblksize)));
             [lam,delta,res] = lanczosfun(Xchol{p},-dX{p},maxit,tol);
             %%
             %% Note: lam <= actual largest eigenvalue <= lam + delta.
             %%
	     d = lam+delta;
          else
             if isempty(invXchol{p});
                invXchol{p} = inv(Xchol{p});
             end
	     tmp = Prod2(pblk,dX{p},invXchol{p},0);
             M = Prod2(pblk,invXchol{p}',tmp,1);
             if (exist('mexblkeig')==3)
                d = mexblkeig(pblk,-M);
             else
                d = blkeig(pblk,-M);
             end
          end
          tmp = max(d) + 1e-15*max(abs(d));
          if (tmp > 0);
             xstep(p) = 1/max(tmp);
          else
             xstep(p) = 1e12;
          end
       elseif strcmp(pblk{1},'q')
          aa = qops(pblk,dX{p},dX{p},2);
          bb = qops(pblk,dX{p},X{p},2);
          cc = qops(pblk,X{p},X{p},2);
          dd = bb.*bb - aa.*cc;
          tmp = min(aa,bb);
          idx = find(dd > 0 & tmp < 0);
          steptmp = 1e12*ones(numblk,1);
          if ~isempty(idx)
             steptmp(idx) = -(bb(idx)+sqrt(dd(idx)))./aa(idx);
          end
          idx = find(abs(aa) < eps & bb < 0);
          if ~isempty(idx)
             steptmp(idx) = -cc(idx)./(2*bb(idx));
          end
          %%
          %% also need first component to be non-negative
          %%
          ss = 1 + [0, cumsum(pblk{2})];
          ss = ss(1:length(pblk{2}));
          dX0 = dX{p}(ss);
          X0 = X{p}(ss);
          idx = find(dX0 < 0 & X0 > 0);
          if ~isempty(idx)
             steptmp(idx) = min(steptmp(idx),-X0(idx)./dX0(idx));
          end
          xstep(p) = min(steptmp);
       elseif strcmp(pblk{1},'l')
          idx = find(dX{p} < 0);
          if ~isempty(idx)
             xstep(p) = min(-X{p}(idx)./dX{p}(idx));
          else
             xstep(p) = 1e12;
          end
       elseif strcmp(pblk{1},'u')
          xstep(p) = 1e12;
       end
    end
    xstep = min(xstep);
%%***************************************************************************
%%***************************************************************************
%% lanczos: find the largest eigenvalue of
%%          invXchol'*dX*invXchol via the lanczos iteration.
%%
%% [lam,delta] = lanczosfun(Xchol,dX,maxit,tol,v)
%%
%% lam:  an estimate of the largest eigenvalue.
%% lam2: an estimate of the second largest eigenvalue.
%% res:  residual norm of the largest eigen-pair.
%% res2: residual norm of the second largest eigen-pair.
%%***************************************************************************

   function [lam,delta,res] = lanczosfun(Xchol,dX,maxit,tol,v)

   if (norm(dX,'fro') < 1e-13)
      lam = 0; delta = 0; res = 0;
      return;
   end
   n = length(dX);
   if (nargin < 5);
      state = randn('state');
      randn('state',0);
      v = randn(n,1);
      randn('state',state);
   end
   if (nargin < 4); tol = 1e-3; end
   if (nargin < 3); maxit = 30; end
   V = zeros(n,maxit+1); H = zeros(maxit+1,maxit);
   v = v/norm(v);
   V(:,1) = v;
   if issparse(Xchol); Xcholtransp = Xchol'; end
%%
%% lanczos iteration.
%%
   for k = 1:maxit
      if issparse(Xchol)
         w = dX*mextriangsp(Xcholtransp,v,1);
         w = mextriangsp(Xchol,w,2);
      else
         w = dX*mextriang(Xchol,v,1);
         w = mextriang(Xchol,w,2);
      end
      wold = w;
      if (k > 1);
         w = w - H(k,k-1)*V(:,k-1);
      end;
      alp = w'*V(:,k);
      w   = w - alp*V(:,k);
      H(k,k) = alp;
      %%
      %% one step of iterative refinement if necessary.
      %%
      if (norm(w) <= 0.8*norm(wold));
         s = (w'*V(:,1:k))';
         w = w - V(:,1:k)*s;
         H(1:k,k) = H(1:k,k) + s;
      end;
      nrm = norm(w);
      v = w/nrm;
      V(:,k+1) = v;
      H(k+1,k) = nrm;  H(k,k+1) = nrm;
      %%
      %% compute ritz pairs and test for convergence
      %%
      if (rem(k,5) == 0) | (k == maxit);
         Hk = H(1:k,1:k); Hk = 0.5*(Hk+Hk');
         [Y,D] = eig(Hk);
         eigH  = real(diag(D));
         [dummy,idx] = sort(eigH);
         res_est = abs(H(k+1,k)*Y(k,idx(k)));
         if (res_est <= 0.1*tol) | (k == maxit);
            lam = eigH(idx(k));
            lam2 = eigH(idx(k-1));
            z  = V(:,1:k)*Y(:,idx(k));
            z2 = V(:,1:k)*Y(:,idx(k-1));
            if issparse(Xchol)
               tmp = dX*mextriangsp(Xcholtransp,z,1);
               res = norm(mextriangsp(Xchol,tmp,2) -lam*z);
               tmp = dX*mextriangsp(Xcholtransp,z2,1);
               res2 = norm(mextriangsp(Xchol,tmp,2) -lam*z2);
            else
               tmp = dX*mextriang(Xchol,z,1);
               res = norm(mextriang(Xchol,tmp,2) -lam*z);
               tmp = dX*mextriang(Xchol,z2,1);
               res2 = norm(mextriang(Xchol,tmp,2) -lam*z2);
            end
            tmp = lam-lam2 -res2;
            if (tmp > 0); beta = tmp; else; beta = eps; end;
            delta = min(res,res^2/beta);
            if (delta <= tol); break; end;
         end
      end
   end
%%***************************************************************************

repo: REST


text:
function y  = laprnd(m, n, mu, sigma)
%LAPRND generate i.i.d. laplacian random number drawn from laplacian distribution
%   with mean mu and standard deviation sigma.
%   mu      : mean
%   sigma   : standard deviation
%   [m, n]  : the dimension of y.
%   Default mu = 0, sigma = 1.
%   For more information, refer to
%   http://en.wikipedia.org./wiki/Laplace_distribution

%   Author  : Elvis Chen (bee33@sjtu.edu.cn)
%   Date    : 01/19/07

%Check inputs
if nargin < 2
    error('At least two inputs are required');
end

if nargin == 2
    mu = 0; sigma = 1;
end

if nargin == 3
    sigma = 1;
end

% Generate Laplacian noise
u = rand(m, n)-0.5;
b = sigma / sqrt(2);
y = mu - b * sign(u).* log(1- 2* abs(u));

author: Sable

repo: mcbench-benchmarks


text:
function c = plus(a,b)
% function C=plus(A,B)
%
% DESCRIPTION
%   Add two polynomials and combine common terms.
%
% INPUTS
%   A,B: polynomials
%
% OUTPUTS
%   C: polynomial, the result of addition A+B
%
% SYNTAX
%   C= A+B
%     Adds A and B.  A and B must have the same dimensions unless one is
%     is a scalar.  A scalar can be added to anything.
%   C = plus(A,B)
%     Function-call form for addition.

% 6/7/2002: PJS  Initial Coding
% 6/8/2002: PJS  Allow matrices of polynomials

% Promote a to polynomial
a = polynomial(a);
sza = size(a);

% Promote b to polynomial
b = polynomial(b);
szb = size(b);

if isempty(a) || isempty(b)

    if isempty(a) && all(szb==[1 1])
        % empty+scalar = empty(sza)
        c=polynomial(zeros(sza));
        return;
    elseif isempty(b) && all(sza==[1 1])
        % scalar+empty = empty(szb)
        c=polynomial(zeros(szb));
        return;
    elseif all(sza==szb)
        c=polynomial(zeros(sza));
        return;
    else
        error('Matrix dimensions must agree.');
    end

elseif all(sza==szb)
    % Matrix + Matrix

    % Get Dimensions
    nta = size(a.degmat,1);
    nva = length(a.varname);
    ntb = size(b.degmat,1);
    nvb = length(b.varname);

    % Stack up coefficients
    acoef = a.coefficient;
    bcoef = b.coefficient;
    coefficient = [acoef; bcoef];

    if nva==0 && nvb==0
        degmat = zeros(nta+ntb,0);
        varname = {};
    else
        adeg = a.degmat;
        bdeg = b.degmat;
        degmat = blkdiag(adeg,bdeg);
        varname = [a.varname(:); b.varname(:)];
    end

    % Form polynomial and combine terms
    chkval = 0; % skip validity check
    c = polynomial(coefficient,degmat,varname,sza,chkval);
    c = combine(c);

elseif all(sza==[1 1]) || all(szb==[1 1])
    % Scalar + Matrix / Matrix + Scalar

    % Make first term of sum be the scalar.
    if ~all(sza==[1 1])
        temp = a;
        a = b;
        b = temp;
        szb = sza;
    end

    % Turn into matrix+matrix
    a.coefficient = repmat(a.coefficient,1,szb(1)*szb(2));
    a.matdim = szb;
    c = plus(a,b);
else
    error('Matrix dimensions must agree');
end

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function result = pyramid_unit_o06_3d ( func )

%*****************************************************************************80
%
%% PYRAMID_UNIT_O06_3D approximates an integral inside the unit pyramid in 3D.
%
%  Discussion:
%
%    A 6 point formula is used.
%
%    The (X,Y,Z) integration region can be represented as:
%
%    - ( 1 - Z ) <= X <= 1 - Z
%    - ( 1 - Z ) <= Y <= 1 - Z
%              0 <= Z <= 1.
%
%    When Z is zero, the integration region is a square lying in the (X,Y)
%    plane, centered at (0,0,0) with "radius" 1.  As Z increases to 1, the
%    radius of the square diminishes, and when Z reaches 1, the square has
%    contracted to the single point (0,0,1).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 April 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carlos Felippa,
%    A compendium of FEM integration formulas for symbolic work,
%    Engineering Computation,
%    Volume 21, Number 8, 2004, pages 867-890.
%
%  Parameters:
%
%    Input, external FUNC, the name of the user supplied function which
%    evaluates F(X,Y,Z), of the form
%      function value = func ( x, y, z )
%
%    Output, real RESULT, the approximate integral of the function.
%
  order = 6;

  w = [ ...
   0.21000000000000000000, ...
   0.21000000000000000000, ...
   0.21000000000000000000, ...
   0.21000000000000000000, ...
   0.06000000000000000000, ...
   0.10000000000000000000 ]';
  x = [ ...
  -0.48795003647426658968, ...
   0.48795003647426658968, ...
   0.48795003647426658968, ...
  -0.48795003647426658968, ...
   0.00000000000000000000, ...
   0.00000000000000000000 ]';
  y = [ ...
  -0.48795003647426658968, ...
  -0.48795003647426658968, ...
   0.48795003647426658968, ...
   0.48795003647426658968, ...
   0.00000000000000000000, ...
   0.00000000000000000000 ]';
  z = [ ...
   0.16666666666666666667, ...
   0.16666666666666666667, ...
   0.16666666666666666667, ...
   0.16666666666666666667, ...
   0.58333333333333333333, ...
   0.75000000000000000000 ]';
%
%  Quadrature.
%
  quad = 0.0;
  for i = 1 : order
    quad = quad + w(i) * feval ( func, x(i), y(i), z(i) );
  end
%
%  Volume.
%
  volume = pyramid_unit_volume_3d ( );
%
%  Result.
%
  result = quad * volume;

  return
end


text:
function classifier = rlsc(options,data)
% {rlsc} trains a Regularized Least Square Classifier (RLSC).
%
%      classifier = rlsc(options,data)
%
%      Note: This function is only a wrapper to 'lapsvmp', where the use of
%            a Hinge function for labeled data is disabled and the
%            contribute of unlabeled data is discarded. See 'lapsvmp'
%            for the description of all required parameters. All parameters
%            related to unlabeled data are obviously not checked.
%
% Author: Stefano Melacci (2009)
%         mela@dii.unisi.it

options.Hinge=0;
options.gamma_I=0;
classifier=lapsvmp(options,data);
classifier.name='rlsc';

author: viggin

repo: domain-adaptation-toolbox


text:
% AddTextToImageWithBorderDemo is a script to demonstrate what
% AddTextToImageWithBorder does.
% It takes about 45 seconds to run on my desktop PC - most of that time is
% spent generating the fonts and displaying the images.
%
% Daniel Warren
% Particle Therapy Cancer Research Institute
% University of Oxford

Image = ones([500 500 3]);

Iterations = 200;

Strings = {'text','image'};

FontNames = {'Arial','Times New Roman','Comic Sans MS','Century Gothic',...
         'Script MT Bold','Courier New','OCR A Extended','Verdana',...
         'Gill Sans MT','Calibri','Cambria'};

FontSizes = 10:20:150;

ThickestBorder = 4;

k = 0;
Fonts = cell(1,length(FontNames)*length(FontSizes));
for i = 1:length(FontNames)
    for j = 1:length(FontSizes)
        k = k+1;
        Fonts{k} = BitmapFont(FontNames{i},FontSizes(j),unique([Strings{:}]));
    end
end

StringIndices = ceil(rand(Iterations,1)*length(Strings));
FontIndices = ceil(rand(Iterations,1)*length(Fonts));
Colors = hsv2rgb([rand(Iterations,1) 0.7+0.3*rand(Iterations,1) 0.9+0.1*rand(Iterations,1)]);
BorderColors = 0.5-(Colors-0.5);
BorderWidths = round(rand(Iterations,1)*ThickestBorder);
XPositions = ceil((-0.2+1.2*rand(1,Iterations))*size(Image,1));
YPositions = ceil((-0.1+1.1*rand(1,Iterations))*size(Image,2));

figure('Units','Pixels','Position',[50 50 50+size(Image,1) 50+size(Image,2)]);
handle = image(Image);
set(gca,'Units','Normalized','Position',[0 0 1 1]);
axis off image;
for i = 1:Iterations
    Image = AddTextToImageWithBorder(Image,Strings{StringIndices(i)},[XPositions(i) YPositions(i)],Colors(i,:),Fonts{FontIndices(i)},[],BorderWidths(i),BorderColors(i,:));
    set(handle,'CData',Image);
    drawnow;
end

author: zouchuhang

repo: LayoutNet


text:
function Z = nonlinearreplace(X,Y,W)

% Very slow, but easily coded. Hopefully not used much

% Special case for replace(x,y,-y)
if isequal(getbase(Y),[0 1]) && isequal(getbase(W),[0 -1]) && isequal(getvariables(Y),getvariables(W)) && is(Y,'linear')
    mt = yalmip('monomtable');
    xv = getvariables(X);
    yv = getvariables(Y);
    if all(mt(:,yv) == fix(mt(:,yv)) )
        signs = ones(length(xv),1);
        for i = 1:length(xv)
            pow = mt(xv(i),yv);
            if ~even(pow)
                signs(i) = -1;
            end
        end
        Z = X;
        Z = setbase(X,getbase(X)*diag([1;signs]));
        return
    end
end

U = sdpvar(length(W),1);
for kk = 1:length(Y)
    Z = [];
    for ii = 1:size(X,1)
        temp = [];
        for jj = 1:size(X,2);

            [coeffs,base] = coefficients(X(ii,jj),Y(kk));

            newp = 0;
            for i = 1:length(base)
                newp = newp + coeffs(i)*U(kk)^degree(base(i));
            end

            temp = [temp newp];
        end
        Z = [Z;temp];
    end
    X = Z;
end

Y = U;
for kk = 1:length(Y)
    Z = [];
    for ii = 1:size(X,1)
        temp = [];
        for jj = 1:size(X,2);

            [coeffs,base] = coefficients(X(ii,jj),Y(kk));

            newp = 0;
            for i = 1:length(base)
                newp = newp + coeffs(i)*W(kk)^degree(base(i));
            end

            temp = [temp newp];
        end
        Z = [Z;temp];
    end
    X = Z;
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [Q, R] = gson(X)
% Gram-Schmidt orthonormalization which produces the same result as [Q,R]=qr(X,0)
% Written by Mo Chen (sth4nth@gmail.com).
[d,n] = size(X);
m = min(d,n);
R = zeros(m,n);
Q = zeros(d,0);
for i = 1:m
    R(1:i-1,i) = Q'*X(:,i);
    v = X(:,i)-Q*R(1:i-1,i);
    R(i,i) = norm(v);
    Q(:,i) = v/R(i,i);
end
R(:,m+1:n) = Q'*X(:,m+1:n);

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
classdef SubcellsMesherFactory < handle

    methods (Access = public, Static)

        function obj = create(cParams)
            ndimIso = cParams.ndimIso;
            type    = cParams.type;
            switch ndimIso
                case 'Line'
                    obj = SubcellsMesher_1D(cParams);
                case 'Surface'
                    switch type
                        case 'INTERIOR'
                            obj  = SubcellsMesher_Interior(cParams);
                        case 'BOUNDARY'
                            obj = SubcellsMesher_Boundary_2D(cParams);
                    end
                case 'Volume'
                    switch type
                        case 'INTERIOR'
                            obj = SubcellsMesher_Interior(cParams);
                        case 'BOUNDARY'
                            obj = SubcellsMesher_Boundary_3D(cParams);
                    end
            end
        end

    end
end


text:
function sR = region(m,varargin)
% return spherical region associated to a set of crystal directions

if check_option(varargin,{'fundamentalRegion','fundamentalSector'}) ...
    && ~check_option(varargin,'complete')
  sR = m.CS.fundamentalSector(varargin{:});
else
  sR = region@vector3d(m,varargin{:});
end

end

author: mtex-toolbox

repo: mtex

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


text:
function [mycat, bigEventCat, bigEventMag] = catalog_overview(in_cat, bigEventMag)
    % catalog_overview presents a window where catalog summary statistics show and can be edited
    % mycat should be a ZmapCatalogView
    %
    %  ask for several input parameters that can be set up
    %  at the beginning of each session. The default values are the
    %  extrema in the catalog
    %
    % bigEventMag typically comes from ZG.CatalogOpts.BigEvents.MinMag

    %ZG=ZmapGlobal.Data;
    report_this_filefun();
    % mycat=ZG.Views.(mycatview);

    %  default values

    %if ~exist('ZG.bin_dur', 'var')   %select bin length respective to time in catalog
    %    ZG.bin_dur = days(30);
    %end
    if isa(in_cat,'ZmapCatalog')
        mycat = ZmapCatalogView(@()in_cat);
    else
        mycat = in_cat;
    end

    big_evt_minmag = ZmapGlobal.Data.CatalogOpts.BigEvents.MinMag;
    bigEventCat = mycat.subset(mycat.Magnitude > big_evt_minmag);
    [minti, maxti] = bounds(mycat.Date);
    daterange = [minti, maxti];
    magrange = bounds2(mycat.Magnitude);
    minma = magrange(1);
    maxma = magrange(2);
    depthrange = bounds2(mycat.Depth);
    mindep = depthrange(1);
    maxdep = depthrange(2);

    fignum = create_dialog();

    watchoff
    %str = 'Please Select a subset of earthquakes and press "Go"';
    %ZmapMessageCenter.set_message('Message',str);
    figure(fignum);

    uiwait(fignum)

    function fignum = create_dialog()
        % create_dialog - creates the dialog box

        %
        % make the interface
        %
        fignum = figure_w_normalized_uicontrolunits(...
            'Units', 'pixels', ...
            'Position',position_in_current_monitor(300, 400),...
            'Name',['Catalog "', mycat.Name,'"'],...
            'Menu','none',...
            'Visible','off',...
            'NumberTitle','off',...
            'Tag', main_dialog_figure('tag'),...
            'NextPlot','new');
        axis off

        % control display parameters
        label_x = 0.05;
        all_h = 0.17;

        %catalog name
        uicontrol('Style','text',...
            'Position',[label_x + 0.01 .94 .3 .05],...
            'Units', 'pixels',...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            ...'FontWeight','bold' ,...
            'HorizontalAlignment','left',...
            'String','Name:');

        uicontrol('Style','edit',...
            'Position',[.4 .94 .5 .05],...
            'Units', 'pixels',...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            ...'FontWeight','bold' ,...
            'HorizontalAlignment','center',...
            'Tag', 'catalog_name_field',...
            'String',mycat.Name,...
            'Callback',@callback_update_catalog_name);

        % EQ's in catalog
        uicontrol('Style','text',...
            'Position',[label_x + 0.01 .88 .7 .05],...
            'Units', 'pixels',...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            ...'FontWeight','bold' ,...
            'HorizontalAlignment','left',...
            'String',' EQs in catalog: ');

        uicontrol('Style','text',...
            'Position',[.70 .88 .22 .05],...
            'Units', 'pixels', ...
            'String',num2str(mycat.Count),...
            'Value',mycat.Count,...
            ...'FontWeight','bold',...
            'FontSize', ZmapGlobal.Data.fontsz.m,...
            'Tag','mapview_nquakes_field',...
            'Callback',@upate_numeric);

        filter_panel = uipanel('Parent',fignum,'Title','Catalog Parameters',...
            'Position',[.02 .4 .96 .45], 'Units', 'pixels');
        option_panel = uipanel('Parent',fignum,'Title','Additional Parameters',...
            'Position',[.02, .15, .96, .25], 'Units', 'pixels');

        % plot big events with M gt
        uicontrol('Parent',option_panel,...
            'Style','text',...
            'Position',[label_x 0.6 .7 .3 ],...
            'Units', 'pixels', ...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            ...'FontWeight','bold' ,...
            'HorizontalAlignment','left',...
            'String','Plot "Big" Events:    M >');

        uicontrol('Parent',option_panel,...
            'Style','edit','Position',[.75 .6 .22 .3],...
            'Units', 'pixels', ...
            'String',num2str(big_evt_minmag),...
            'Value',big_evt_minmag,...
            'Tag','mapview_big_evt_field',...
            'Callback',@update_editfield_value);

        % TODO: add reset button (would be nice...)

        %  beginning year
        uicontrol('Parent',filter_panel,...
            'Style','edit','Position',[.45 0.75 .52 all_h],...
            'Units','pixels',...
            'String',char(minti,'uuuu-MM-dd HH:mm:ss'),...
            'Value', datenum(minti),...
            'HorizontalAlignment','center',...
            'Callback',@update_editfield_date,...
            'Tag','mapview_start_field',...
            'tooltipstring', 'as decimal year or yyyy-mm-dd hh:mm:ss');

        uicontrol('Parent',filter_panel,...
            'Style','text',...
            'Position',[label_x 0.75 .395 all_h],...
            'Units', 'pixels', ...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            ...'FontWeight','bold' ,...
            'HorizontalAlignment','left',...
            'String','Beginning date: ');

        % ending year
        uicontrol('Parent',filter_panel,...
            'Style','edit','Position',[.45 0.55 .52 all_h],...
            'Units', 'pixels', ...
            'String',char(maxti,'uuuu-MM-dd HH:mm:ss'),...
            'Value', datenum(maxti),...
            'Callback',@update_editfield_date,...
            'Tag','mapview_end_field',...
            'tooltipstring', 'as decimal year or yyyy-mm-dd hh:mm:ss');

        uicontrol('Parent',filter_panel,...
            'Style','text',...
            'Position',[label_x 0.55 0.4 all_h],...
            'Units', 'pixels', ...
            'FontSize',ZmapGlobal.Data.fontsz.m ,...
            ...'FontWeight','bold' ,...
            'HorizontalAlignment','left',...
            'String','Ending date: ');

        % Magnitude

        uicontrol('Parent',filter_panel,...
            'Style','text',...
            'Position',[label_x 0.29 0.5 all_h],...
            'Units', 'pixels', ...
            'FontSize',ZmapGlobal.Data.fontsz.m,...
            ...'FontWeight','bold',...
            'HorizontalAlignment','left',...
            'String','Magnitude:');

        uicontrol('Parent',filter_panel,...
            'Style','text',...
            'Position',[0.63 0.29 0.4 all_h],...
            'Units', 'pixels', ...
            'FontSize',ZmapGlobal.Data.fontsz.m,...
            ...'FontWeight','bold',...
            'HorizontalAlignment','left',...
            'String','<= M <=');

        uicontrol('Parent',filter_panel,...
            'Style','edit',...
            'Position',[.43 .3 .17 all_h],...
            'Units', 'pixels', ...
            'String',num2str(minma),...
            'Value', minma,...
            'Tag','mapview_minmag_field',...
            'Callback',@update_editfield_value);

        uicontrol('Parent',filter_panel,...
            'Style','edit','Position',[.80 .3 .17 all_h],...
            'Units', 'pixels', ...
            'String',num2str(maxma),...
            'Value', maxma,...
            'Tag','mapview_maxmag_field',...
            'Callback',@update_editfield_value);

        % Depth control
        uicontrol('Parent',filter_panel,...
            'Style','text',...
            'Position',[label_x 0.04 0.4 all_h],...
            'Units', 'pixels', ...
            'FontSize',ZmapGlobal.Data.fontsz.m,...
            ...'FontWeight','bold',...
            'HorizontalAlignment','left',...
            'String','Depth:');

        uicontrol('Parent',filter_panel,...
            'Style','text',...
            'Position',[0.5 0.04 0.4 all_h],...
            'Units', 'pixels', ...
            'FontSize',ZmapGlobal.Data.fontsz.m,...
            ...'FontWeight','bold',...
            'HorizontalAlignment','left',...
            'String','<= Z (km) <=');

        uicontrol('Parent',filter_panel,...
            'Style','edit','Position',[.29 .05 .2 all_h],...
            'Units', 'pixels', ...
            'String',num2str(mindep),...
            'Value', mindep,...
            'Tag','mapview_mindepth_field',...
            'Callback',@update_editfield_value);

        uicontrol('Parent',filter_panel,...
            'Style','edit','Position',[.77 .05 .2 all_h],...
            'Units', 'pixels', ...
            'String',num2str(maxdep),...
            'Value', maxdep,...
            'Tag','mapview_maxdepth_field',...
            'Callback',@update_editfield_value);

        % buttons
        uicontrol('Style','Pushbutton',...
            'Position',[.79 .02 .20 .10 ],...
            'Units', 'pixels', ...
            'Callback',@cancel_callback,'String','Cancel');

        uicontrol('Style','Pushbutton',...
            'Position',[.58 .02 .20 .10 ],...
            'Units','pixels',...
            'Callback',@go_callback,...
            'String','Apply');

        uicontrol('Style','Pushbutton',...
            'Position',[.05 .02 .35 .10 ],...
            'Units', 'pixels', ...
            'Callback',{@distro_callback, mycat}, ...
            'String','see distributions');

        set(gcf,'visible','on')
    end

    function go_callback(src, ~)
        %TODO remove all the side-effects.  store relevent data somewhere specific
        %filter the catalog, then return
        disp(mycat)
        ZG = ZmapGlobal.Data;
        myparent=get(src,'Parent');

        h = findall(myparent,'Tag','mapview_maxdepth_field');
        maxdep = h.Value;
        h = findall(myparent,'Tag','mapview_minmag_field');
        minma = h.Value;
        h = findall(myparent,'Tag','mapview_maxmag_field');
        maxma = h.Value;
        h = findall(myparent,'Tag','mapview_mindepth_field');
        mindep = h.Value;
        h = findall(myparent,'Tag','mapview_start_field');
        minti = datetime(datevec(h.Value));
        h = findall(myparent,'Tag','mapview_end_field');
        maxti = datetime(datevec(h.Value));
        h = findall(myparent,'Tag','mapview_big_evt_field');
        ZG.CatalogOpts.BigEvents.MinMag = h.Value;
        h = findall(myparent,'Tag','catalog_name_field');
        mycat.Name = h.String;
        %h = findall(myparent,'Tag','mapview_binlen_field');
        %ZG.bin_dur = days(h.Value);

        % following code originally from sele_sub.m
        %    Create  reduced (in time and magnitude) catalogs "primeCatalog" and "newcat"
        %
        mycat.DepthLims=[mindep, maxdep];
        mycat.DateLims=[minti, maxti];
        mycat.MagnitudeLims=[minma, maxma];

        %create catalog of "big events" if not merged with the original one:
        %
        bigEventCat = mycat.subset(mycat.Magnitude > ZG.CatalogOpts.BigEvents.MinMag);

        %mycat.sort('Date');
        % ZG.Views.(mycatview)=mycat;

        % ZmapMessageCenter.update_catalog();
        %zmap_update_displays();

        close(main_dialog_figure('handle'));
        disp('changed to...')
        disp(mycat);
        % changes in bin length go to global
    end

    function callback_update_catalog_name(~,~)
        disp('in name callback')
        % if field should respond somehow to changes, do it here
    end

end

function cancel_callback(~, ~)
    % return without making changes to catalog
    ZmapMessageCenter.update_catalog();
    %h=ZmapMessageCenter();
    %h.update_catalog();
    close(main_dialog_figure('handle'));
end

function distro_callback(src,~,mycat)
    watchon; drawnow;
    dlg = main_dialog_figure('handle');
    if numel(dlg) >1
        warning('ZMAP:nonUniqueHandle','multiple dialog windows found')
    end

    f = findall(0,'Tag','catoverview_distribution_pane');
    if isempty(f)
        % grow the catalog figure. create plots in the empty portion, change the button behavior
        dlg.Position = dlg.Position + [0 0 450 0];
        src.String = 'hide distributions';
        pp=uipanel('Parent',gcf,'Units','pixels','Position',[310 10 420 370],'Tag','catoverview_distribution_pane','Title','Distributions');
        %f = figure('Name','Catalog time, mag, depth distributions','MenuBar','none','NumberTitle','off','Tag','catoverview_distribution');
        t_p=subplot(3,1,1,'Parent',pp);
        m_p=subplot(3,1,2,'Parent',pp);
        d_p=subplot(3,1,3,'Parent',pp);
        histogram(t_p,mycat.Date);
        xlabel(t_p,'Time');
        histogram(m_p,mycat.Magnitude);
        xlabel(m_p,'Magnitude');
        histogram(d_p,mycat.Depth);
        xlabel(d_p,'Depth');
    else
        delete(findobj('Tag','catoverview_distribution_pane'));
        dlg.Position = dlg.Position - [0 0 450 0];

        src.String = 'show distributions';
        % delete the histograms,
    end
    watchoff; drawnow;
end

%% Tag Name Helpers
function answer = main_dialog_figure(opt)
    % get
    % opt is either 'tag' or 'handle'
    s = 'catalog_overview_dlg';
    switch opt
        case 'tag'
            answer = s;
        case 'handle'
            answer = findobj('Tag', s);
        otherwise
            error("main_dialog_figure:invalid option, must be 'tag' or 'handle'");
    end
end

author: CelsoReyes


text:
function [B,D,mu] = extract_params_from_gbn(bnet)
% Extract all the local parameters of each Gaussian node, and collect them into global matrices.
% [B,D,mu] = extract_params_from_gbn(bnet)
%
% B(i,j) is a block matrix that contains the transposed weight matrix from node i to node j.
% D(i,i) is a block matrix that contains the noise covariance matrix for node i.
% mu(i) is a block vector that contains the shifted noise mean for node i.

% In Shachter's model, the mean of each node in the global gaussian is
% the same as the node's local unconditional mean.
% In Alag's model (which we use), the global mean gets shifted.

num_nodes = length(bnet.dag);
bs = bnet.node_sizes(:); % bs = block sizes
N = sum(bs); % num scalar nodes

B = zeros(N,N);
D = zeros(N,N);
mu = zeros(N,1);

for i=1:num_nodes % in topological order
  ps = parents(bnet.dag, i);
  e = bnet.equiv_class(i);
  %[m, Sigma, weights] = extract_params_from_CPD(bnet.CPD{e});
  s = struct(bnet.CPD{e}); % violate privacy of object
  m = s.mean; Sigma = s.cov; weights = s.weights;
  if length(ps) == 0
    mu(block(i,bs)) = m;
  else
    mu(block(i,bs)) = m + weights *  mu(block(ps,bs));
  end
  B(block(ps,bs), block(i,bs)) = weights';
  D(block(i,bs), block(i,bs)) = Sigma;
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function kern = tensorKernParamInit(kern)

% TENSORKERNPARAMINIT TENSOR kernel parameter initialisation.
% The tensor product (TENSOR) kernel is a container kernel for
% allowing tensor products kernels of separate component kernels.
% several different kernels to be added together. It is created by
% using the kernCreate command with the kernel type given as a
% cell and the first entry given as 'tensor'. For example, to
% create a tensor kernel compound kernel that is composed of
% an RBF kernel and a LIN kernel you call,
%
% kern = kernCreate(X, {'tensor', 'rbf', 'lin'});
%
% Each individual kernel is then stored within the returned kernel
% structure. The kernels are stored in order in a field called
% 'comp'. So display obtain the 'rbf' kernel you write:
%
% kernDisplay(kern.comp{1})
%
% SEEALSO : cmpndKernParamInit, multiKernParamInit
%
% FORMAT
% DESC initialises the tensor product
%  kernel structure with some default parameters.
% ARG kern : the kernel structure which requires initialisation.
% RETURN kern : the kernel structure with the default parameters placed in.
%
% SEEALSO : kernCreate, kernParamInit
%
% COPYRIGHT : Neil D. Lawrence, 2006

% KERN

kern.nParams = 0;
kern.transforms = [];
if ~isfield(kern, 'comp')
  kern.comp=cell(0);
end
for i = 1:length(kern.comp)
  kern.comp{i} = kernParamInit(kern.comp{i});
  kern.nParams = kern.nParams + kern.comp{i}.nParams;
  kern.comp{i}.index = [];
end
kern.paramGroups = speye(kern.nParams);

% Warn if component kernels have white variance and check if
% resulting kernel is stationary.
kern.isStationary = true;
whiteVariance = 0;
for i = 1:length(kern.comp)
  if ~kern.comp{i}.isStationary
    kern.isStationary = false;
  end
  if strcmp(kern.comp{i}.type, 'white')
    whiteVariance = whiteVariance + kern.comp{i}.variance;
  else
    if(isfield(kern.comp{i}, 'whiteVariance'))
      whiteVariance = whiteVariance + ...
          kern.comp{i}.whiteVariance;
    end
  end
end
if whiteVariance > 0
  warning(['Components of tensor kernel have non-zero white variance. ' ...
            'This can cause problems some implementations.']);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [A, C, b, f] = HALS(Y, A, C, b, f, params)
%% Hierarchical alternating least square method for solving NMF problem
% Y = A*C + b*f

%input:
%   Y:      d1 X d2 X T, raw data. It will be reshaped to (d1*d2) X T in this
%   function
%   A:      (d1*d2) X K, initial value of spatial components
%   C:      K X T, initial value of temporal components
%   b:      (d1*d2) X nb, initial value of background spatial component
%   f:      nb X T, initial value of background temporal component
%   params: parameters used in this function.
%       bSiz:   blur size. A box kernel (bSiz X bSiz) will be convolved
%       with each neuron's initial spatial component, then all nonzero
%       pixels will be picked as pixels to be updated, and the rest will be
%       forced to be 0.
%       maxIter: maximum iteration of iterating HALS.

% Author: Pengcheng Zhou, Carnegie Mellon University, based on a python
% implementation from Johannes Friedrich, Columbia University, 2015.

%% parameters
if isfield(params, 'maxIter'), maxIter = params.maxIter; else maxIter=5; end
if isfield(params, 'search_method'); method=params.search_method; else method='ellipse'; end
if and(isfield(params, 'bSiz'), strcmpi(method, 'dilate'))
    params.se = strel('disk', params.bSiz);
end
Y = reshape(Y, size(A, 1), []);
% search locations
IND = determine_search_location(A, method, params);

%% update spatial and temporal components neuron by neurons

for miter=1:maxIter
    %% update neurons
    Yac = Y - b*f;
    ind_del = find(std(A,0,1)==0);
    A(:, ind_del) = [];
    C(ind_del, :) = [];
    IND(:, ind_del) = [];
    %   temporal
    C = HALS_temporal(Yac, A, C, 5);

    ind_del = find(std(C,0,2)==0);
    A(:, ind_del) = [];
    C(ind_del, :) = [];
    IND(:, ind_del) = [];
    %   spatial
    A = HALS_spatial(Yac, A, C, IND, 5);

    %% update background
    Ybg = Y-A*C;
    % temporal
    f = HALS_temporal(Ybg, b, f, 5);
    % spatial
    b = HALS_spatial(Ybg, b, f, [], 5);
end

author: zhoupc


text:
% HISTWC  Weighted histogram count given number of bins
%
% This function generates a vector of cumulative weights for data
% histogram. Equal number of bins will be considered using minimum and
% maximum values of the data. Weights will be summed in the given bin.
%
% Usage: [histw, vinterval] = histwc(vv, ww, nbins)
%
% Arguments:
%       vv    - values as a vector
%       ww    - weights as a vector
%       nbins - number of bins
%
% Returns:
%       histw     - weighted histogram
%       vinterval - intervals used
%
%
%
% See also: HISTC, HISTWCV

% Author:
% mehmet.suzen physics org
% BSD License
% July 2013

function [histw, vinterval] = histwc(vv, ww, nbins)
  minV  = min(vv);
  maxV  = max(vv);
  delta = (maxV-minV)/nbins;
  vinterval = linspace(minV, maxV, nbins)-delta/2.0;
  histw = zeros(nbins, 1);
  for i=1:length(vv)
    ind = find(vinterval < vv(i), 1, 'last' );
    if ~isempty(ind)
      histw(ind) = histw(ind) + ww(i);
    end
  end
end


text:
%% AMG TEST IV: Three Dimensional Problem
%
% We consider the linear finite element discretization of Poisson equation
% in three dimensions with Dirichlet boundary conditions and compare
% geometric multigrid and algebraic multigrid.

%% Conclusion
%
% The algebraic multigrid is robust to the size of the matrix. The
% iteration steps increases slightly. The preformance is better on
% structure grids than unstructure grids. The sparsity of the coarse grid
% is increased.

%% Uniform Mesh: Geometric Multigrid
clear all
option.solver = 'mg';
cubePoisson;

%% Uniform Mesh: Algebraic Multigrid
clear all
option.solver = 'amg';
cubePoisson;

%% Unstructured Mesh
clear all;
load oilpump;
showboundary3(node,elem);
option.solver = 'amg';
cubePoisson;

repo: Aero-matlab


text:
function [ x, w ] = nco_nested_rule ( n )

%*****************************************************************************80
%
%% NCO_NESTED_RULE computes the nested Newton Cotes Open rule.
%
%  Discussion:
%
%    Nested version of Newton Cotes Open rule on [-1,+1].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 February 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the rule.
%
%    Output, real X(N), W(N), the points and weights.
%

%
%  Get N equally spaced points that do not include [ 0, 1];
%
  x = vdc_sequence ( n );
%
%  Linear transformation from [0,1] to [-1,+1].
%
  x = 2.0 * x - 1.0;
%
%  Compute the weights.
%
  x_min = -1.0;
  x_max = +1.0;

  w = nc_compute ( n, x_min, x_max, x );

  return
end


text:
function [ n_data, nu, x, fx ] = bessel_ix_values ( n_data )

%*****************************************************************************80
%
%% BESSEL_IX_VALUES returns some values of the Ix Bessel function.
%
%  Discussion:
%
%    This set of data considers the less common case in which the
%    index of the Bessel function In is actually not an integer.
%    We may suggest this case by occasionally replacing the symbol
%    "In" by "Ix".
%
%    The modified Bessel functions In(Z) and Kn(Z) are solutions of
%    the differential equation
%
%      Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0.
%
%    In Mathematica, the function can be evaluated by:
%
%      BesselI[n,x]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Cambridge University Press, 1999,
%    ISBN: 0-521-64314-7,
%    LC: QA76.95.W65.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real NU, the order of the function.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 28;

  fx_vec = [ ...
    0.3592084175833614E+00,  ...
    0.9376748882454876E+00,  ...
    2.046236863089055E+00,   ...
    3.053093538196718E+00,   ...
    4.614822903407601E+00,   ...
    26.47754749755907E+00,   ...
    2778.784603874571E+00,   ...
    4.327974627242893E+07,   ...
    0.2935253263474798E+00,  ...
    1.099473188633110E+00,   ...
    21.18444226479414E+00,   ...
    2500.906154942118E+00,   ...
    2.866653715931464E+20,   ...
    0.05709890920304825E+00, ...
    0.3970270801393905E+00,  ...
    13.76688213868258E+00,   ...
    2028.512757391936E+00,   ...
    2.753157630035402E+20,   ...
    0.4139416015642352E+00,  ...
    1.340196758982897E+00,   ...
    22.85715510364670E+00,   ...
    2593.006763432002E+00,   ...
    2.886630075077766E+20,   ...
    0.03590910483251082E+00, ...
    0.2931108636266483E+00,  ...
    11.99397010023068E+00,   ...
    1894.575731562383E+00,   ...
    2.716911375760483E+20 ];

  nu_vec = [ ...
    0.50E+00, ...
    0.50E+00, ...
    0.50E+00, ...
    0.50E+00, ...
    0.50E+00, ...
    0.50E+00, ...
    0.50E+00, ...
    0.50E+00, ...
    1.50E+00, ...
    1.50E+00, ...
    1.50E+00, ...
    1.50E+00, ...
    1.50E+00, ...
    2.50E+00, ...
    2.50E+00, ...
    2.50E+00, ...
    2.50E+00, ...
    2.50E+00, ...
    1.25E+00, ...
    1.25E+00, ...
    1.25E+00, ...
    1.25E+00, ...
    1.25E+00, ...
    2.75E+00, ...
    2.75E+00, ...
    2.75E+00, ...
    2.75E+00, ...
    2.75E+00 ];

  x_vec = [ ...
      0.2E+00, ...
      1.0E+00, ...
      2.0E+00, ...
      2.5E+00, ...
      3.0E+00, ...
      5.0E+00, ...
     10.0E+00, ...
     20.0E+00, ...
      1.0E+00, ...
      2.0E+00, ...
      5.0E+00, ...
     10.0E+00, ...
     50.0E+00, ...
      1.0E+00, ...
      2.0E+00, ...
      5.0E+00, ...
     10.0E+00, ...
     50.0E+00, ...
      1.0E+00, ...
      2.0E+00, ...
      5.0E+00, ...
     10.0E+00, ...
     50.0E+00, ...
      1.0E+00, ...
      2.0E+00, ...
      5.0E+00, ...
     10.0E+00, ...
     50.0E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    nu = 0.0;
    x = 0.0;
    fx = 0.0;
  else
    nu = nu_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function datapt = lvd_StopwatchMassTasks(stateLogEntry, subTask, stopwatch)
%lvd_StopwatchMassTasks Summary of this function goes here
%   Detailed explanation goes here

    switch subTask
        case 'swValue'
            stopwatchStates = stateLogEntry.getAllStopwatchStates();
            stopwatchState = stopwatchStates([stopwatchStates.stopwatch] == stopwatch);

            if(not(isempty(stopwatchState)))
                stopwatchState = stopwatchState(1);
                datapt = stopwatchState.value;
            else
                datapt = -1;
            end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function h = spm_orthviews_showposition()
% :Usage:
% ::
%
%    handles = spm_orthviews_showposition;
%
% plots x, y, and z coordinates on SPM orthviews figure
%
% bring context structure variable into this script.
% contains handles, etc.
%
% ..
%    tor wager, august 2006
%    updated: june 2007
% ..

    global st

    % get current position in mm
    pos = round(spm_orthviews('Pos'));  % x y z

    for i = 1:length(st.vols)

        if ~isempty(st.vols{i}) && ~isempty(st.vols{i}.ax) && ~isempty(st.vols{i}.ax{1})

            axish(1) = st.vols{i}.ax{1}.ax;
            axish(2) = st.vols{i}.ax{2}.ax;
            axish(3) = st.vols{i}.ax{3}.ax;

            % axial
            h(1) = put_text(axish(1),['z = ' num2str(pos(3))]);

            % coronal
            h(2) = put_text(axish(2),['y = ' num2str(pos(2))]);

            % saggital
            h(3) = put_text(axish(3),['x = ' num2str(pos(1))]);

        end

        % %
        % % % Get figure and axis handles
        % % fh = findobj('Tag','Graphics');
        % % ch = get(fh,'Children');
        % % for i= 1:length(ch),mytype = get(ch(i),'Type'); wh(i)=strcmp(mytype,'axes'); end
        % %
        % % if ~exist('wh','var')
        % %     disp('spm_orthviews_showposition: No axis handles found. Exiting.');
        % %     return
        % % end
        % %
        % % axish = ch(find(wh));
        % %
        % % if isempty(axish), return, end

        % % % get which axis is which
        % % for i = 1:length(axish)
        % % poss(i,:) = get(axish(i),'Position');
        % % end
        % %
        % % % get rid of extra axes we may have created in the 4th quadrant
        % % other_axes = find(any(poss(:,1) > .45 & poss(:,2) < .2,2));
        % % axish(other_axes) = [];
        % % poss(other_axes,:) = [];
        % %
        % % % sort into order:  axial, coronal, saggital
        % % ssum = sum(poss(:,1:2),2);
        % % [ssum,ind] = sort(ssum);
        % % axish = axish(ind);
        % %
        % %
        % % % axial
        % % h(1) = put_text(axish(1),['z = ' num2str(pos(3))]);
        % %
        % % % coronal
        % % h(2) = put_text(axish(2),['y = ' num2str(pos(2))]);
        % %
        % % % saggital
        % % h(3) = put_text(axish(3),['x = ' num2str(pos(1))]);

    end  % end loop thru vols

    return

function h = put_text(axish,str)

    axes(axish)

    % remove old text
    oldh = findobj(gca,'Type','text');
    delete(oldh)

    xdelta = diff(get(axish,'XLim'));
    ydelta = diff(get(axish,'YLim'));

    % add new text
    xpos = xdelta .* .42;
    ypos = ydelta .* .1;
    h1 = text(xpos - xdelta * .01,ypos + ydelta * .005,str,'Color','k','FontSize',24, 'FontWeight', 'b');
    h = text(xpos,ypos,str,'Color','g','FontSize',24, 'FontWeight', 'demi');

    return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%% Specifications
minspec=70;
maxspec=200;

%% Find
found= data(:)>minspec & data(:)<maxspec;
Percent=sum(found)/length(found)


text:
function stackOut = ComputeLDRStackHistogram(stack)
%
%       stackOut = ComputeLDRStackHistogram(stack)
%
%
%        Input:
%           -dir_name: the folder name where the stack is stored as a
%           series of LDR images.
%           -format: an LDR format for reading LDR images in the current directory
%
%        Output:
%           -stackOut: a stack of LDR image histograms
%
%     Copyright (C) 2013  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

[~, ~, col, n] = size(stack);

stackOut = zeros(256, col, n);

for i=1:n
    %store in the stack
    for j=1:col

        tmp = stack(:,:,j,i);

        if(isa(tmp, 'double') || isa(tmp, 'single'))
            tmp = uint8(ClampImg(round(tmp * 255), 0.0, 255.0));
        end

        if(isa(tmp, 'uint16'))
            tmp = uint8(ClampImg(round(tmp / 255), 0, 255));
            warning('Is this a 16-bit image? The maximum is set to 65535.');
        end

        stackOut(:,j,i) = imhist(tmp);
    end
end

end

author: banterle


text:
function idx = findClosestCentroids(X, centroids)
%FINDCLOSESTCENTROIDS computes the centroid memberships for every example
%   idx = FINDCLOSESTCENTROIDS (X, centroids) returns the closest centroids
%   in idx for a dataset X where each row is a single example. idx = m x 1
%   vector of centroid assignments (i.e. each entry in range [1..K])
%

% Set K
K = size(centroids, 1);

% You need to return the following variables correctly.
idx = zeros(size(X,1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every example, find its closest centroid, and store
%               the index inside idx at the appropriate location.
%               Concretely, idx(i) should contain the index of the centroid
%               closest to example i. Hence, it should be a value in the
%               range 1..K
%
% Note: You can use a for-loop over the examples to compute this.
%
for i = 1:size(X,1)
	min = Inf;
	for j = 1:K
		diff = sum((X(i,:) - centroids(j,:)).^2);
		if min > diff
			min = diff;
			idx(i) = j;
		end
	end
end

% =============================================================

end

author: AvaisP


text:
function t = dec2twos(x, nbits)

% DEC2TWOS Convert decimal integer to binary string two's complement.
%
% Usage: T = DEC2TWOS(X, NBITS)
%
% Converts the signed decimal integer given by X (either a scalar, vector,
% or matrix) to the two's complement representation as a string. If X is a
% vector or matrix, T(i, :) is the representation for X(i) (the shape of X
% is not preserved).
%
% Example:
%     >> dec2twos([23 3 -23 -3])
%
%     ans =
%
%     010111
%     000011
%     101001
%     111101
%
% Inputs:
%   -X: decimal integers to convert to two's complement.
%   -NBITS: number of bits in the representation (optional, default is the
%   fewest number of bits necessary).
%
% Outputs:
%   -T: two's complement representation of X as a string.
%
% See also: TWOS2DEC, DEC2FIX, FIX2DEC, DEC2BIN, DEC2HEX, DEC2BASE.

error(nargchk(1, 2, nargin));
x = x(:);
maxx = max(abs(x));
nbits_min = nextpow2(maxx + (any(x == maxx))) + 1;

% Default number of bits.
if nargin == 1 || isempty(nbits)
    nbits = nbits_min;
elseif nbits < nbits_min
    warning('dec2twos:nbitsTooSmall', ['Minimum number of bits to ' ...
        'represent maximum input x is %i, which is greater than ' ...
        'input nbits = %i. Setting nbits = %i.'], ...
        nbits_min, nbits, nbits_min)
    nbits = nbits_min;
end

t = repmat('0', numel(x), nbits); % Initialize output:  Case for x = 0
if any(x > 0)
    t(x > 0, :) = dec2bin(x(x > 0), nbits);           % Case for x > 0
end
if any(x < 0)
    t(x < 0, :) = dec2bin(2^nbits + x(x < 0), nbits); % Case for x < 0
end


text:
%QUAD2TRI Convert quadrilaterals to triangles.
%
%   [ SOUT ] = QUAD2TRI( SIN, IORDER ) Convert quadrilateral to triangular cells.
%
%       Input       Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       sin         struct                 Grid or problem struct with
%                                          p, c (and optionally s and b) fields
%       iorder      1/2/{0}                Node ordering of triangles.
%                                            0 : Automatically select best splitting
%                                            1 : 1-2-3 / 1-3-4 splitting
%                                            2 : 2-3-4 / 1-2-4 splitting
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       sout        struct                 Output grid or problem struct
%
%   See also HEX2TET, TET2HEX, TRI2QUAD

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function op = linop_subsample( sz, omega, SYMMETRIC )
%LINOP_SUBSAMPLE Subsampling linear operator.
%OP = LINOP_SUBSAMPLE( SZ, OMEGA )
%   vector and matrix subsampling. Depending on SZ and OMEGA,
%   this can do row-sampling (e.g. a partial FFT)
%   or it can sample specific entries of a matrix (e.g. matrix completion)
%   To specify matrix subsampling, set sz = { [n1,n2], [length(omega),1] }
%       where the input matrix has size n1 x n2
%
%OP = LINOP_SUBSAMPLE( OMEGA )
%   works if OMEGA is a sparse matrix, whose nonzero entries
%   specify the entries to sample.
%
%OP = LINOP_SUBSAMPLE( ..., SYMMETRIC )
%   If SYMMETRIC is true, then
%   forces the domain to be the space of symmetric matrices
%
% Example with FFT:
%   FFT = linop_handles([N,N], @(x)fft(x)/sqrt(N), @(x)sqrt(N)*real(ifft(x)) ,'R2C');
%   and to make a row sub-sampled FFT (M <= N):
%   rows = randperm(M);
%   A   = linop_compose( linop_subsample([M,N], rows ), FFT );

% Designed to be used with a fft or dct
% Do we need a separate oversampling version?
% Also, make linop_DCT and linop_FFT ?
%   The reason we might want this is that for linop_FFT,
%   people will probably use idct as the transpose -- this is not
%   correct, due to scaling

% Help documentation: TBD
%    Constructs a TFOCS-compatible linear operator from separate function
%    handles that compute the forward and adjoint operations. The first
%    argument, SZ, gives the size of the linear operator; and the forward
%    and adjoint handles are AF and AT, respectively.
%
%    If the inputs and outputs are simple vectors, then SZ can take the
%    standard Matlab form [N,M], where N is the input length and M is the
%    output length. If the input or output is a matrix or array, then SZ
%    should take the form { S_in, S_out }, where S_in is the size of the
%    input and S_out is the size of the output.
%

error(nargchk(1,3,nargin));
if nargin < 3
    SYMMETRIC = false;
end
if nargin == 1
    omega = sz;
    sz = { size(omega), [nnz(omega),1] };
elseif nargin ==2 && issparse(sz)
    SYMMETRIC = omega;
    omega = sz;
    sz = { size(omega), [nnz(omega),1] };
end

if numel( sz ) ~= 2,
    error( 'Size must have two elements.' );
elseif isnumeric( sz ),
    sz = { [sz(2),1], [sz(1),1] };
elseif ~isa( sz, 'cell' ),
    error( 'Invalid operator size specification.' );
end

dom = sz{1};
n1 = dom(1); n2 = dom(2);
ran = sz{2};

% There are several possibilities.  Let x be a point in the domain
%{
    x is a vector.  Then omega should be a vector.
                    We return x(omega), resp.

    x is a matrix.
        omega is a vector
            This is ambiguous: does the user want x(omega,:) or x(omega)?
        omega is a matrix with 2 columns, then assume it is [I,J]
            We convert it to linear indices.
        omega is a general matrix, more than 2 columns
            Not sure what the user means; report an error.
        omega is a sparse matrix
            We find it's entries, and use those

%}

if n2 == 1
    % x is a vector, not a matrix.  Simple.
%     op = @(x,mode) linop_subsample_vector( sz, omega, x, mode );
    % Allow it to vectorize along rows:
    op = @(x,mode) linop_subsample_row( sz, omega, x, mode );
else
    % trickier case.
    if issparse(omega)
        ind = find(omega);
        [I,J] = ind2sub( sz{1}, ind );
%         op = @(x,mode) linop_subsample_matrix( sz, ind, I, J, SYMMETRIC, x, mode );
    elseif isvector(omega)
        ind     = omega;
        [I,J]   = ind2sub( sz{1}, ind );
%         op = @(x,mode) linop_subsample_matrix( sz, omega, SYMMETRIC,x, mode );
    elseif size(omega,2) == 2
        ind = sub2ind( sz{1}, omega(:,1), omega(:,2) );
        I   = omega(:,1);
        J   = omega(:,2);
%         op = @(x,mode) linop_subsample_matrix( sz, ind, SYMMETRIC,x, mode );
    else
        error('omega is not an acceptable size; perhaps you meant it to be sparse?');
    end
    % Make sure 'ind' is a column vector
    ind     = ind(:);
    op = @(x,mode) linop_subsample_matrix( sz, ind, I, J, SYMMETRIC, x, mode );

end

function y = linop_subsample_vector(sz, omega, x, mode )
switch mode,
    case 0, y = sz;
    case 1,
        if ~isequal( sz{1}, size(x) ), error('input wrong size for vector subsampling'); end
        y = x(omega);
    case 2,
        y = zeros( sz{1} );
        y(omega) = x;
end
function y = linop_subsample_row(sz, omega, x, mode )
switch mode,
    case 0, y = sz;
    case 1,
        if ~isequal( sz{1}(1), size(x,1) ), error('input wrong size for vector subsampling'); end
        y = x(omega,:);
    case 2,
        y = zeros( sz{1} );
        y = zeros( sz{1}(1), size(x,2) );
        y(omega,:) = x;
end
function y = linop_subsample_matrix(sz, omega, indI, indJ,SYMMETRIC, x, mode )
switch mode,
    case 0, y = sz;
    case 1
        S = [];
        S.type = '()';
        S.subs = {omega};
        y = subsref(x,S);
    case 2,
        dom = sz{1}; n1 = dom(1); n2 = dom(2);
        y = sparse( indI, indJ, x, n1, n2 );
        if SYMMETRIC
            % in future, might update this
            % e.g. force omega to only refer to lower part of matrx,
            % and then do the update y = y + tril(y,-1)'
            y = (y+y')/2;
        end
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function [ r, z, c, s ] = schex ( r, ldr, p, k, l, z, ldz, nz, job )

%*****************************************************************************80
%
%% SCHEX updates the Cholesky factorization of a positive definite matrix.
%
%  Discussion:
%
%    The factorization has the form
%
%      A = R' * R
%
%    where A is a positive definite matrix of order P.
%
%    The updating involves diagonal permutations of the form
%
%      E' * A * E
%
%    where E is a permutation matrix.  Specifically, given
%    an upper triangular matrix R and a permutation matrix
%    E (which is specified by K, L, and JOB), SCHEX determines
%    an orthogonal matrix U such that
%
%      U * R * E = RR,
%
%    where RR is upper triangular.  At the user's option, the
%    transformation U will be multiplied into the array Z.
%    If A = X'*X, so that R is the triangular part of the
%    QR factorization of X, then RR is the triangular part of the
%    QR factorization of X*E, that is, X with its columns permuted.
%
%    For a less terse description of what SCHEX does and how
%    it may be applied, see the LINPACK guide.
%
%    The matrix Q is determined as the product U(L-K)*...*U(1)
%    of plane rotations of the form
%
%      (    C(I)       S(I) )
%      (                    ),
%      (   -S(I)       C(I) )
%
%    where C(I) is real, the rows these rotations operate on
%    are described below.
%
%    There are two types of permutations, which are determined
%    by the value of JOB.
%
%    1, right circular shift.  The columns are rearranged in the order:
%
%         1,...,K-1,L,K,K+1,...,L-1,L+1,...,P.
%
%       U is the product of L-K rotations U(I), where U(I)
%       acts in the (L-I,L-I+1)-plane.
%
%    2, left circular shift: the columns are rearranged in the order
%
%         1,...,K-1,K+1,K+2,...,L,K,L+1,...,P.
%
%       U is the product of L-K rotations U(I), where U(I)
%       acts in the (K+I-1,K+I)-plane.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real R(LDR,P), the upper triangular factor that is to be updated.
%    Elements of R below the diagonal are not referenced.
%
%    Input, integer LDR, the leading dimension of the array R.
%    LDR must be at least P.
%
%    Input, integer P, the order of the matrix R.
%
%    Input, integer K, the first column to be permuted.
%
%    Input, integer L, the last column to be permuted.
%    L must be strictly greater than K.
%
%    Input, real Z(LDZ,NZ), an array of NZ P-vectors into
%    which the transformation U is multiplied.  Z is not referenced if NZ = 0.
%
%    Input, integer LDZ, the leading dimension of the array Z.
%    LDZ must be at least P.
%
%    Input, integer NZ, the number of columns of the matrix Z.
%
%    Input, integer JOB, determines the type of permutation.
%    1, right circular shift.
%    2, left circular shift.
%
%    Output, real R(LDR,P), the updated upper triangular factor.
%
%    Output real Z(LDZ,NZ), the updated array of vectors.
%
%    Output, real C(P), S(P), the cosines and sines of the
%    transforming rotations.
%

%
%  Initialize
%
  lmk = l - k;
  lm1 = l - 1;
%
%  Right circular shift.
%
  if ( job == 1 )
%
%  Reorder the columns.
%
    for i = 1 : l
      ii = l - i + 1;
      s(i) = r(ii,l);
    end

    for jj = k : lm1
      j = lm1 - jj + k;
      for i = 1 : j
        r(i,j+1) = r(i,j);
      end
      r(j+1,j+1) = 0.0;
    end

    for i = 1 : k-1
      ii = l - i + 1;
      r(i,k) = s(ii);
    end
%
%  Calculate the rotations.
%
    t = s(1);
    for i = 1 : lmk
      [ c(i), s(i), s(i+1), t ] = srotg ( s(i+1), t );
      t = s(i+1);
    end

    r(k,k) = t;

    for j = k+1 : p
      il = max ( 1, l-j+1 );
      for ii = il : lmk
        i = l - ii;
        t = c(ii) * r(i,j) + s(ii) * r(i+1,j);
        r(i+1,j) = c(ii) * r(i+1,j) - s(ii) * r(i,j);
        r(i,j) = t;
      end
    end
%
%  If required, apply the transformations to Z.
%
    for j = 1 : nz
      for ii = 1 : lmk
        i = l - ii;
        t = c(ii) * z(i,j) + s(ii) * z(i+1,j);
        z(i+1,j) = c(ii) * z(i+1,j) - s(ii) * z(i,j);
        z(i,j) = t;
      end
    end
%
%  Left circular shift.
%
  else
%
%  Reorder the columns.
%
    for i = 1 : k
      ii = lmk + i;
      s(ii) = r(i,k);
    end

    for j = k : lm1
      for i = 1 : j;
        r(i,j) = r(i,j+1);
      end
      jj = j - k + 1;
      s(jj) = r(j+1,j+1);
    end

    for i = 1 : k
      ii = lmk + i;
      r(i,l) = s(ii);
    end

    r(k+1:l,l) = 0.0;
%
%  Reduction loop.
%
    for j = k : p
%
%  Apply the rotations.
%
      if ( j ~= k )

        iu = min ( j-1, l-1 );

        for i = k : iu
          ii = i - k + 1;
          t = c(ii) * r(i,j) + s(ii) * r(i+1,j);
          r(i+1,j) = c(ii) * r(i+1,j) - s(ii) * r(i,j);
          r(i,j) = t;
        end

      end

      if ( j < l )
        jj = j - k + 1;
        t = s(jj);
        [ c(jj), s(jj), r(j,j), t ] = srotg ( r(j,j), t );
      end

    end
%
%  Apply the rotations to Z.
%
    for j = 1 : nz
      for i = k : lm1
        ii = i - k + 1;
        t = c(ii) * z(i,j) + s(ii) * z(i+1,j);
        z(i+1,j) = c(ii) * z(i+1,j) - s(ii) * z(i,j);
        z(i,j) = t;
      end
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m