text:
function [kappa,pval,k, pk]=kappa_stats(D,ncat)
% calculate Kappa statisc for agreement
% D: matrix NxM of N scores by M raters
% cat: categories
% Written By Issam El Naqa Date: 03/21/07
% ref.: SAS docs. (Fleiss '81)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

[N,M]=size(D);
lk=length(ncat); % number of categories
for i=1:lk
    x(:,i) = sum(D==ncat(i),2);
end
p=sum(x)/(N*M);
k=1-sum(x.*(M-x))./((N*M*(M-1).*p.*(1-p))+eps);
sek=sqrt(2/(N*M*(M-1)));
pk=drxlr_get_p_gaussian(k./sek)/2; % one-sided
% kappa=sum(p.*(1-p).*k)/sum(p.*(1-p));
kappa=1-(N*M^2-sum(sum(x.^2)))/((N*M*(M-1)*sum(p.*(1-p)))+eps);
sekappa=sqrt(2)/(sum(p.*(1-p)*sqrt(N*M*(M-1)))+eps)*sqrt(sum(p.*(1-p))^2-sum(p.*(1-p).*(1-2*p)));
z=kappa/sekappa;
pval=drxlr_get_p_gaussian(z)/2; % one-sided
return

function p=drxlr_get_p_gaussian(x)
% two tailed p-value from a normal distribution
%DREX subfunction
%Written by Issam El Naqa 2003-2005
%Extracted for generalized use 2005, AJH

    p = erfc(abs(x)./sqrt(2));

return

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function result = polygon_1_2d ( n, v )

%*****************************************************************************80
%
%% POLYGON_1_2D integrates the function 1 over a polygon in 2D.
%
%  Discussion
%
%    The polygon is bounded by the points (X(1:N), Y(1:N)).
%
%    INTEGRAL = 0.5 * SUM ( 1 <= I <= N )
%      ( X(I) + X(I-1) ) * ( Y(I) - Y(I-1) )
%
%    where X(0) and Y(0) should be replaced by X(N) and Y(N).
%
%    Note that the integral of 1 over a polygon is the area of the polygon.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    S F Bockman,
%    Generalizing the Formula for Areas of Polygons to Moments,
%    American Mathematical Society Monthly,
%    1989, pages 131-132.
%
%  Parameters:
%
%    Input, integer N, the number of vertices of the polygon.
%    N should be at least 3 for a nonzero result.
%
%    Input, real V(2,N), the coordinates of the vertices
%    of the polygon.  These vertices should be given in counter-clockwise order.
%
%    Output, real RESULT, the value of the integral.
%
  result = 0.0;

  if ( n < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'POLYGON_1_2D - Warning!\n' );
    fprintf ( 1, '  The number of vertices must be at least 3.\n' );
    fprintf ( 1, '  The input value of N = %d\n', n );
    error ( 'POLYGON_1_2D - Fatal error!' );
  end

  for i = 1 : n

    if ( i == 1 )
      im1 = n;
    else
      im1 = i - 1;
    end

    result = result + 0.5 * ( v(1,i) + v(1,im1) ) * ( v(2,i) - v(2,im1) );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function chebyshev_polynomial_test16 ( )

%*****************************************************************************80
%
%% CHEBYSHEV_POLYNOMIAL_TEST16 tests W_POLYNOMIAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHEBYSHEV_POLYNOMIAL_TEST16:\n' );
  fprintf ( 1, '  W_POLYNOMIAL_VALUES stores values of\n' );
  fprintf ( 1, '  the Chebyshev polynomials.\n' );
  fprintf ( 1, '  W_POLYNOMIAL evaluates the polynomial.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '                        Tabulated                 Computed\n' );
  fprintf ( 1, '     N        X           W(n,x)                    W(n,x)                     Error\n' );

  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, x, fx1 ] = w_polynomial_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2_vec = w_polynomial ( 1, n, x );
    fx2 = fx2_vec(1,n+1);
    e = fx1 - fx2;

    fprintf ( 1, '  %4d  %12f  %24.16e  %24.16e  %8.2g\n', n, x, fx1, fx2, e );

  end

  return
end


text:
function x = binary(x)
%BINARY Constrains variables to be binary (0/1)
%
%   F = BINARY(x) is used to a posteriori constrain
%   variables to be binary, in contrast to BINVAR
%   that declares variables as binary a priori.
%
%   NOTE
%    The binary constraint is imposed on the involved
%    decision variables, not on the actual SDPVAR object.
%
%   INPUT
%    x : SDPVAR object
%
%   OUTPUT
%    F : Constraint object
%
%   EXAMPLE
%    F = binary(x);                  % Add integrality
%    F = binary(x*pi)                % Equivalent to code above
%
%   See also INTEGER, SET, SDPVAR, INTVAR, BINVAR

x.typeflag = 8;
x = lmi(x);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function bvec2 = bvec_reverse ( n, bvec1 )

%*****************************************************************************80
%
%% BVEC_REVERSE reverses a binary vector.
%
%  Discussion:
%
%    A BVEC is an integer vector of binary digits, intended to
%    represent an integer.  BVEC(1) is the units digit, BVEC(N-1)
%    is the coefficient of 2**(N-2), and BVEC(N) contains sign
%    information.  It is 0 if the number is positive, and 1 if
%    the number is negative.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 November 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the length of the vectors.
%
%    Input, integer BVEC1(N), the vector to be reversed.
%
%    Output, integer BVEC2(N), the reversed vector.
%
  bvec2(1:n) = bvec1(n:1:-1);

  return
end


text:
function [features, labels, labelprobs, weights] = ...
    mcmcGetAllGoodSegments(imsegs, imdir, vclassifierSP, hclassifierSP, ...
    eclassifier, segclassifier, spdata, adjlist, edata, ncv)

for f = 1:numel(imsegs)
    k = mod(ceil(f/(numel(imsegs)/ncv)), ncv) + 1;
    disp([num2str(f) ': ' imsegs(f).imname])
    im = im2double(imread([imdir '/' imsegs(f).imname]));
    [features{f}, labels{f}, labelprobs{f}, weights{f}] = ...
        mcmcGenerateGoodSegments(im, imsegs(f), vclassifierSP, hclassifierSP, ...
        eclassifier, segclassifier(k), spdata{f}, adjlist{f}, edata{f});
end

author: zouchuhang

repo: LayoutNet


text:
function [bler, ber] = simulation(N, K, design_epsilon, max_runs, max_err, resolution, ebno_vec, list_size_vec, g, crc_length)
R = (K - crc_length)/N;
lambda_offset = 2.^(0 : log2(N));
llr_layer_vec = get_llr_layer(N);
bit_layer_vec = get_bit_layer(N);
[G_crc, H_crc] = crc_generator_matrix(g, K - crc_length);
crc_parity_check = G_crc';

%beta expansion Code Construction, proposed by Huawei
% beta = sqrt(sqrt(2));
% % channels = PW(N, beta);
% % channels = HPW(N, beta);
% channels = EPW(N, beta);
% [~, channel_ordered] = sort(channels, 'descend');
% info_bits = sort(channel_ordered(1 : K));
% frozen_bits = ones(N , 1);
% frozen_bits(info_bits) = 0;

%Bhattacharyya Code (BEC) Construction
% channels = get_BEC_IWi(N, design_epsilon);
% [~, channel_ordered] = sort(channels, 'descend');
% info_bits = sort(channel_ordered(1 : K), 'ascend');
% frozen_bits = ones(N , 1);
% frozen_bits(info_bits) = 0;
% info_bits_logical = logical(mod(frozen_bits + 1, 2));

%Gaussian approximation Code Construction
design_snr = 2.5;
sigma_cc = 1/sqrt(2 * R) * 10^(-design_snr/20);
[channels, ~] = GA(sigma_cc, N);
[~, channel_ordered] = sort(channels, 'descend');
info_bits = sort(channel_ordered(1 : K), 'ascend');
frozen_bits = ones(N, 1);
frozen_bits(info_bits) = 0;
frozen_bits = logical(frozen_bits);

%Following MATLAB code is used for matrix-multiplication based systematic polar encoding
FN = get_GN(N);
GAA = FN(info_bits, info_bits);
GAAC = FN(info_bits, frozen_bits);
Generate_xAC = mod(GAAC' * GAA', 2);
systematic_encoding_algorithm = 1;

%This value can be 1, 2, 3, and 4.
%1 : 'Matrix_multiplication_systematic_encoder'
%!!!!But when N is large,the Generator matrix will be large, too. The computer may not have such large storage.
%2 : 'Sarkis_Two_step_systematic_encoder'
%3 : 'Arikan_SC_style_systematic_encoder'
%4 : 'Arikan_Recursive_systematic_encoder'. This one is slow.

%Special constituent nodes
node_type_matrix = get_node_structure(frozen_bits);
psi_vec = get_psi_for_advanced_sc_decoder(node_type_matrix);

%Results Stored
bler = zeros(length(ebno_vec), length(list_size_vec));
num_runs = zeros(length(ebno_vec), length(list_size_vec));
ber = zeros(length(ebno_vec), length(list_size_vec));
%Loop starts
tic
for i_run = 1 : max_runs
    if  mod(i_run, max_runs/resolution) == 1
        disp(' ');
        disp(['Sim iteration running = ' num2str(i_run)]);
        disp(['N = ' num2str(N) ' K = ' num2str(K)]);
        disp(['List size = ' num2str(list_size_vec)]);
        disp('Current block error performance');
        disp(num2str([ebno_vec'  bler./num_runs]));
        disp('Current bit error performance');
        disp(num2str([ebno_vec'  ber./num_runs/K]));
        disp(' ')
    end
    info  = rand(K - crc_length, 1) > 0.5;
    info_with_crc = [info; mod(crc_parity_check * info, 2)];
    switch systematic_encoding_algorithm
        case 1
            parity_check_bits = mod(Generate_xAC * info_with_crc, 2);
            x = zeros(N, 1);
            x(info_bits) = info_with_crc;
            x(frozen_bits) = parity_check_bits;
        case 2
            x = sarkis_systematic_polar_encoder(info_with_crc, info_bits, frozen_bits, N, lambda_offset, llr_layer_vec);
        case 3
            x = arikan_sc_systematic_polar_encoder(info_with_crc, frozen_bits, info_bits, lambda_offset, llr_layer_vec, bit_layer_vec);
        case 4
            x = zeros(N, 1);
            x(info_bits) = info_with_crc;
            x = arikan_recursive_systematic_polar_encoder(x, mod(frozen_bits + 1, 2));
    end
    bpsk = 1 - 2 * x;
    noise = randn(N, 1);
    prev_decoded = zeros(length(ebno_vec), length(list_size_vec));
    for i_ebno = 1 : length(ebno_vec)
        sigma = 1/sqrt(2 * R) * 10^(-ebno_vec(i_ebno)/20);
        y = bpsk + sigma * noise;
        llr = 2/sigma^2*y;
        %*******Simulaion Accelaration*********
        for i_list = 1 : length(list_size_vec)
            if i_list ~= 1
                if bler(i_ebno, i_list) == max_err
                    continue;
                end
            else
                if all(bler(i_ebno, 2 : end) == max_err)
                    continue
                end
            end
            num_runs(i_ebno, i_list) = num_runs(i_ebno, i_list) + 1;
            run_sim = 1;
            for i_ebno2 = 1 : i_ebno
                for i_list2 = 1 : i_list
                    if prev_decoded(i_ebno2, i_list2)
                        run_sim = 0;
                    end
                end
            end
            if run_sim == 0
                continue;
            end
            %*******Simulaion Accelaration*********
            if list_size_vec(i_list) == 1
                polar_info_esti = SC_decoder(llr, frozen_bits, lambda_offset, llr_layer_vec, bit_layer_vec, info_bits);
%                 polar_info_esti = FastSCdecoder(llr, info_bits, node_type_matrix, lambda_offset, llr_layer_vec, psi_vec, bit_layer_vec);
            else
                polar_info_esti = FastSCL_decoder(llr, list_size_vec(i_list), info_bits, lambda_offset, llr_layer_vec, bit_layer_vec, psi_vec, node_type_matrix, H_crc);
            end
            if any(polar_info_esti ~= info_with_crc)
                bler(i_ebno, i_list) = bler(i_ebno, i_list) + 1;
                ber(i_ebno, i_list) = ber(i_ebno, i_list) + sum(polar_info_esti(1 : K - crc_length) ~= info_with_crc(1 : K - crc_length));
                %Caution, this is data bits error rate, regardless of CRC
                %bits.
            else
                prev_decoded(i_ebno, i_list) = 1;
            end

        end
    end
end
toc
end

author: YuYongRun

repo: PolarCodeDecodersInMatlab


text:
function [rtk,ins_align_flag]=ins_align(rtk,obsr_,obsb_,nav)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2020-2025, by Kai Chen, All rights reserved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global glc
ins_align_flag=0;

if norm(rtk.sol.pos)~=0
    [vel,ins_align_flag]=tdcp2vel(rtk,nav,obsr_,rtk.oldobsr);
end

[rtk,~]=gnss_solver(rtk,obsr_,obsb_,nav);

rtk.oldobsr=obsr_;

if rtk.sol.stat~=glc.SOLQ_NONE&&ins_align_flag==1
    % initialize ins
    pos=xyz2blh(rtk.sol.pos);
    yaw=vel2yaw(vel);
    att=[0 0 yaw];
    avp0=[att,vel,pos]';
    ins=ins_init(rtk.opt.ins,avp0);

    % correct lever arm for position and velocity
    ins.pos=ins.pos-ins.Mpv*ins.Cnb*ins.lever;
    ins.vel=ins.vel-ins.Cnb*askew(ins.web)*ins.lever;
    rtk.ins=ins;
else
    ins_align_flag=0;
end

return

repo: GINav


text:
function test_homer

% WALLTIME 00:10:00
% MEM 2GB
% DEPENDENCY homer2fieldtrip fieldtrip2homer ft_write_data opto2homer homer2opto

% there are different ways that Homer data can be processed with FieldTrip
% - the way that is used in most FT tutorials is to read (original) Homer files using FT_PREPROCESSING
% - an alternative way is to read (original) Homer files using homer2fieldtrip
% - it is also possible to write FieldTrip data to a Homer file on disk

filename = dccnpath('/home/common/matlab/fieldtrip/data/test/original/nirs/homer/S1001_run01.nirs');

data1 = homer2fieldtrip(filename);
data2 = homer2fieldtrip(fieldtrip2homer(data1));

p = tempdir;
f1 = fullfile(p, 'data1.nirs');
f2 = fullfile(p, 'data2.nirs');

ft_write_data(f1, data1.trial{1}, 'header', data1.hdr)
ft_write_data(f2, data2.trial{1}, 'header', data2.hdr)

% the files written by FieldTrip should still have the same content
data1a = homer2fieldtrip(f1);
data2a = homer2fieldtrip(f2);

% reading the files using the tutorial approach with FT_DEFINETRIAL and FT_PREPROCESSING
cfg = [];
cfg.dataset = filename;
data0b = ft_preprocessing(cfg);
cfg.dataset = f1;
data1b = ft_preprocessing(cfg);
cfg.dataset = f2;
data2b = ft_preprocessing(cfg);

% clean up
delete(f1)
delete(f2)

%%

% do not compare the hdr field, it is different in hdr.orig
fn = {'time', 'trial', 'label', 'opto'};

for i=1:numel(fn)
  field = fn{i};
  % after roundtrip fieldtrip->homer->fieldtrip
  assert(isalmostequal(data1.(field), data2. (field), 'abstol', 1e-8), sprintf('%s is different between data1 and data2', field));

  % after ft_write_data
  assert(isalmostequal(data1.(field), data1a.(field), 'abstol', 1e-8), sprintf('%s is different between data1 and data1a', field));
  assert(isalmostequal(data1.(field), data2a.(field), 'abstol', 1e-8), sprintf('%s is different between data1 and data2a', field));

  % after ft_preprocessing
  assert(isalmostequal(data1.(field), data0b.(field), 'abstol', 1e-8), sprintf('%s is different between data1 and data0b', field));
  assert(isalmostequal(data1.(field), data1b.(field), 'abstol', 1e-8), sprintf('%s is different between data1 and data1b', field));
  assert(isalmostequal(data1.(field), data2b.(field), 'abstol', 1e-8), sprintf('%s is different between data1 and data2b', field));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% This subroutine assigns creates a grid with
% spacing dx,dy (in degreees). The size will
% be selected interactiVELY. The Dvalue in each
% volume around a grid point containing ni earthquakes
% will be calculated. This code is called from dcparain.m.
%
%   Stefan Wiemer 1/95

global no1 bo1 inb1 inb2 eq0p

% the new data vector to be analysed is called Da, relative to the conter of the x-section and already in km
% D = [x,y,z ]
Da = [eq0p(1,:)' eq0p(2,:)' a.Date a.Date.Month a.Date.Day a.Magnitude a.Depth];
Da0 = find(Da(:,7) > 0);
Da = Da.subset(Da0);
clear Da0;

if sel == 'in'
    % get the grid parameter
    % initial values
    %
    dd = 1.00;
    dx = 1.00 ;
    ni = 600;
    Nmin = 600; %on line 303 it has been replaced by ni
    stan2 = nan;
    stan = nan;
    prf = nan;
    av = nan;

    % make the interface
    %
    figure_w_normalized_uicontrolunits(...
        'Name','Grid Parameters',...
        'NumberTitle','off', ...
        'MenuBar','none', ...
        'units','points',...
        'Visible','on', ...
        'Position',[ 100 200 500 200]);
    axis off
    % Francesco ...

    % creates a dialog box to input grid parameters
    %
    freq_field=uicontrol('Style','edit',...
        'Position',[.32 .57 .12 .08],...
        'Units','normalized','String',num2str(ni),...
        'Callback','ni=str2double(get(freq_field,''String'')); set(freq_field,''String'',num2str(ni));set(tgl2,''value'',0); set(tgl1,''value'',1)');

    freq_field0=uicontrol('Style','edit',...
        'Position',[.80 .57 .12 .08],...
        'Units','normalized','String',num2str(ra),...
        'Callback','ra=str2double(get(freq_field0,''String'')); set(freq_field0,''String'',num2str(ra)) ; set(tgl2,''value'',1); set(tgl1,''value'',0)');

    freq_field2=uicontrol('Style','edit',...
        'Position',[.32 .44 .12 .08],...
        'Units','normalized','String',num2str(dx),...
        'Callback','dx=str2double(get(freq_field2,''String'')); set(freq_field2,''String'',num2str(dx));');

    freq_field3=uicontrol('Style','edit',...
        'Position',[.32 .31 .12 .08],...
        'Units','normalized','String',num2str(dd),...
        'Callback','dd=str2double(get(freq_field3,''String'')); set(freq_field3,''String'',num2str(dd));');

    tgl1 = uicontrol('Backgroundcolor', [0.8 0.8 0.8], 'Fontweight','bold',...
        'FontSize', 10, 'Style','checkbox',...
        'string','Number of Events:',...
        'Position',[.05 .56 .2 .10],...
         'Callback','set(tgl2,''value'',0, ''ForegroundColor'', ''w''); set(tgl1, ''ForegroundColor'', ''k'')',...
        'Units','normalized');

    set(tgl1,'value',1);

    tgl2 =  uicontrol('BackGroundColor', [0.8 0.8 0.8],'Style','checkbox',...
        'string','Constant Radius:','Fontweight','bold','FontSize', 10,...
        'Position',[.55 .56 .2 .1],...
         'Callback','set(tgl1,''value'',0,''ForegroundColor'', ''w''); set(tgl2, ''ForegroundColor'', ''k'')',...
        'Units','normalized');

    set(tgl2, 'ForegroundColor', 'w');

    close_button=uicontrol('Style','Pushbutton',...
        'Position',[.45 .05 .15 .13 ],...
        'Units','normalized','Callback','close;done','String','Cancel');

    help_button=uicontrol('Style','Pushbutton',...
        'Position',[.70 .05 .15 .13 ],...
        'Units','normalized','Callback','close;done','String','Help');

    go_button1=uicontrol('Style','Pushbutton',...
        'Position',[.20 .05 .15 .13 ],...
        'Units','normalized',...
        'Callback','tgl1 =get(tgl1,''Value'');tgl2 =get(tgl2,''Value'');close; sel = ''ca''; rdwstest; ',...
        'String','Go');

    txt3 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0.35 0.9 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.l ,...
        'FontWeight','bold',...
        'String',' Grid Parameters');
    txt5 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[-0.07 0.46 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Horizontal Spacing [km]:');

    txt6 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[-0.07 0.30 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Depth spacing [km]:');

    txt7 = text(...
        'Color',[0 0 0 ],...
        'EraseMode','normal',...
        'Position',[0.45 0.62 0 ],...
        'Rotation',0 ,...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'Color', 'r',...
        'String','OR');

    set(gcf,'visible','on');
    watchoff

end   % if sel == in

% get the grid-size interactively and
% calculate the b-value in the grid by sorting
% the seimicity and selectiong the ni neighbors
% to each grid point

if sel == 'ca'

    figure_w_normalized_uicontrolunits(xsec_fig)
    hold on

    messtext=...
        ['To select a polygon for a grid.       '
        'Please use the LEFT mouse button of   '
        'or the cursor to the select the poly- '
        'gon. Use the RIGTH mouse button for   '
        'the final point.                      '
        'Mac Users: Use the keyboard "p" more  '
        'point to select, "l" last point.      '
        '                                      '];

    zmap_message_center.set_message('Select Polygon for a grid',messtext);

    figure;
    ax=plot(Da(:,1),-Da(:,7),'o');
    xlabel('Distance in [km]')
    ylabel('Depth in [km]')

    hold on
    %ax = findobj('Tag','main_map_ax');
    [x,y, mouse_points_overlay] = select_polygon(ax);
    zmap_message_center.set_info('Message',' Thank you .... ')

    plos2 = plot(x,y,'b-','era','xor', 'Color', 'r');        % plot outline
    sum3 = 0.;
    pause(0.3)

    %create a rectangular grid
    xvect=[min(x):dx:max(x)];
    yvect=[min(y):dd:max(y)];
    gx = xvect;gy = yvect;
    tmpgri=zeros((length(xvect)*length(yvect)),2);
    n=0;

    for i=1:length(xvect)
        for j=1:length(yvect)
            n=n+1;
            tmpgri(n,:)=[xvect(i) yvect(j)];
        end
    end

    %extract all gridpoints in chosen polygon
    XI=tmpgri(:,1);
    YI=tmpgri(:,2);

    ll = polygon_filter(x,y, XI, YI, 'inside');
    %grid points in polygon
    newgri=tmpgri(ll,:);

    % Plot all grid points
    plot(newgri(:,1),newgri(:,2),'+k','era','back')

    if length(xvect) < 2  ||  length(yvect) < 2
        errordlg('Selection too small! (not a matrix)');
        return
    end

    itotal = length(newgri(:,1));
    if length(gx) < 4  ||  length(gy) < 4
        errordlg('Selection too small! ');
        return
    end

    zmap_message_center.set_info(' ','Running... ');think
    %  make grid, calculate start- endtime etc.  ...
    %
    t0b = min(newa.Date)  ;
    n = newa.Count;
    teb = newa(n,3) ;
    tdiff = round((teb - t0b)*365/par1);
    loc = zeros(3, length(gx)*length(gy));

    % loop over  all points
    %
    i2 = 0.;
    i1 = 0.;
    wai = waitbar(0,' Please Wait ...  ');
    set(wai,'NumberTitle','off','Name','D-value grid - percent done', 'position', [250 80 270 50]);
    drawnow;
    %
    %
    %
    sdrd = [];
    k = 1;
    for ni = 50:50:1000

        rdtest = [];

        for i= 1:length(newgri(:,1))
            x = newgri(i,1);y = newgri(i,2);
            %allcount = allcount + 1.;
            i2 = i2+1;

            % calculate distance from center point and sort wrt distance
            % Francesco: please check if the distance  l is right

            l = sqrt(((Da(:,1) - x)).^2 + ((Da(:,7) + y)).^2 + (Da(:,2).^2)) ;
            [s,is] = sort(l);
            b = Da(is(:,1),:) ;       % re-orders matrix to agree row-wise

            if tgl1 == 0   % take point within r
                l3 = l <= ra;
                b = Da.subset(l3);      % new data per grid point (b) is sorted in distanc
                rd = ra;
            else
                % take first ni points
                b = b(1:ni,:);      % new data per grid point (b) is sorted in distance
                l2 = sort(l);
                rd = l2(ni);

            end

            %estimate the completeness and b-value, and take the zero depth events away.
            %newt2 = [b(:,1) b(:,2) zeros(size(b,1),1) zeros(size(b,1),1) zeros(size(b,1),1) zeros(size(b,1),1) b(:,3)];
            %
            newt2  = b;
            E = newt2;

            if length(b) >= ni  % enough events?

                rdtest = [rdtest; rd];

            end %if length >= ni

        end  % for  newgri

        mrd = sum(rdtest,2)/size(rdtest,2);
        vrd = var(rdtest);
        %varrd1 = (rdtest-mrd).^2;
        %varrd = sum(varrd1, 1)*1/size(rdtest,1);
        stdevrd = vrd.^(1/2);
        sdrd = [sdrd; stdevrd];

        waitbar((1/20)*k);
        k = k+1;
    end %for ni

    clear mrd vrd
    close(wai)

    sdrd
    radius1 = find(min(sdrd));
    radius = (radius1*50) + 50

    clear radius1

end %if sel = ca

author: CelsoReyes


text:
function sws = sw_normalized(sw, swc, sor)
%SW_NORMALIZED maps the saturation from [0, 1] to [swc, 1-sor]
sws=((sw>swc).*(sw<1-sor).*(sw-swc)./(1-sor-swc)+(sw>=1-sor));
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function As = slsymgraph(A, symmethod)
%SLSYMGRAPH Forces symmetry of the adjacency matrix of a graph
%
% $ Syntax $
%   - As = slsymgraph(A)
%   - As = slsymgraph(A, symmethod)
%
% $ Arguments $
%   - A:            The adjacency matrix of the original graph
%   - symmethod:    The method to symmetrize the graph
%   - As:           The symmetry adjacency matrix
%
% $ Description $
%   - As = slsymgraph(A) makes a symmetry version of the adjacency matrix
%     A using default method.
%
%   - As = slsymgraph(A, symmethod) symmetrizes the adjacency matrix by
%     using the specified method.
%     \*
%     \t    Table. The method to symmetrizes adjacency matrix
%     \h       name       &              description
%             'avgor'     & Force symmetry using the following rule:
%                           if both aij and aji are non-zeros, then
%                           take their average
%                           if only one of aij and aji is non-zero, then
%                           take the non-zero one
%                           if both aij and aji are zeros, then set zero
%                           (for both logical and numeric)
%             'avgand'    & Force symmetry using the following rule:
%                           if both aij and aji are non-zeros, then
%                           take their average
%                           if either one of aij and aji is zero, then
%                           set zero
%                           (for both logical and numeric)
%             'or'        & Use or-rule: d = aij | aji
%                           (for only logical)
%             'and'       & Use and-rule: d = aij & aji
%                           (for only logical)
%             'simavg'    & make simple average: always take (aij+aji)/2
%                           (for only numeric)
%     \*
%     The default method to use is 'avgor'. You can use your own function
%     handle. It should be like the following form:
%       v = f(v1, v2)
%     v1 and v2 are arrays of equal size with corresponding values. For
%     a reasonable function, it should satisfy that:
%       f(v, v) == v && f(v1, v2) = f(v2, v1)
%
% $ Remarks $
%   - A can be full matrix or sparse matrix, and As preserves the same
%     storage form.
%
%   - A should be a square matrix.
%
%   - When 'avgor' method is applied to logical, it is equivalent to 'or'
%     when 'avgand' method is applied to logical, it is equivalent to
%     'and'.
%
% $ History $
%   - Created by Dahua Lin, on Sep 8, 2006
%

%% parse and verify input

if ndims(A) ~= 2 || size(A,1) ~= size(A,2)
    error('sltoolbox:invalidarg', ...
        'The A should be a square 2D matrix');
end
n = size(A, 1);

if nargin < 2 || isempty(symmethod)
    symmethod = 'avgor';
end

if ischar(symmethod)
    switch symmethod
        case 'avgor'
            fcs = @compsym_avgor;
        case 'avgand'
            fcs = @compsym_avgand;
        case 'or'
            fcs = @compsym_or;
            if isnumeric(A)
                error('sltoolbox:rterror', ...
                    'The or method is not applicable to numerical matrix');
            end
        case 'and'
            fcs = @compsym_and;
            if isnumeric(A)
                error('sltoolbox:rterror', ...
                    'The and method is not applicable to numerical matrix');
            end
        otherwise
            error('sltoolbox:invalidarg', ...
                'Invalid method for symmetrization: %s', method);
    end
elseif isa(symmethod, 'function_handle')
    fcs = symmethod;
else
    error('sltoolbox:invalidarg', ...
        'Invalid method for symmetrization.');
end

%% main skeleton

% prepare all indices with aij or aji non-zero

[I0, J0] = find(A);

% single out diagonal ones
is_diag = (I0 == J0);
if any(is_diag)
    inds_diag = sub2ind([n, n], I0(is_diag), J0(is_diag));
else
    inds_diag = [];
end

% process the non-diagonal ones
not_diag = ~is_diag;
clear is_diag;

if any(not_diag)
    % filter indices
    I0 = I0(not_diag);
    J0 = J0(not_diag);
    clear not_diag;

    % merge to down triangular part
    I = I0;
    J = J0;
    idx_ut = find(I0 > J0);
    if ~isempty(idx_ut)
        I(idx_ut) = J0(idx_ut);
        J(idx_ut) = I0(idx_ut);
    end
    clear I0 J0 idx_ut;

    % unique and expand to up triangular part
    inds_dt = sub2ind([n, n], I, J);
    [inds_dt, si] = unique(inds_dt);
    I = I(si);
    J = J(si);
    inds_ut = sub2ind([n, n], J, I);
    clear I J si;
else
    inds_dt = [];
    inds_ut = [];
end

% get original values

if ~isempty(inds_dt)
    v_dt = A(inds_dt);
    v_ut = A(inds_ut);
else
    v_dt = [];
    v_ut = [];
end

% compute the symmetrized value

v = fcs(v_dt, v_ut);
clear v_dt v_ut;

% combine diagonal value and non-diagonal value

if ~isempty(inds_diag)
    v_diag = A(inds_diag);
else
    v_diag = [];
end

s_inds = vertcat(inds_diag, inds_dt, inds_ut);
clear inds_diag inds_dt inds_ut;
s_vals = vertcat(v_diag, v, v);
clear v_diag v;

% create matrix

As = slmakeadjmat(n, n, s_inds, s_vals, islogical(A), issparse(A));

%% symmetry value computation functions

function vd = compsym_avgor(v1, v2)

if isnumeric(v1)
    has_both = v1 & v2;
    only_v1 = v1 & ~v2;
    only_v2 = v2 & ~v1;

    vd = zeros(size(v1));
    vd(has_both) = (v1(has_both) + v2(has_both)) / 2;
    vd(only_v1) = v1(only_v1);
    vd(only_v2) = v2(only_v2);
else
    vd = v1 | v2;
end

function vd = compsym_avgand(v1, v2)

if isnumeric(v1)
    has_both = v1 & v2;

    vd = zeros(size(v1));
    vd(has_both) = (v1(has_both) + v2(has_both)) / 2;
else
    vd = v1 & v2;
end

function vd = compsym_or(v1, v2)

vd = v1 | v2;

function vd = compsym_and(v1, v2)

vd = v1 & v2;

author: lmthang


text:
function xfm = xfmPoints(x,y)

% add third dimension
x = [x zeros(size(x,1),1)];
x(1,3) =1;
y = [y ones(size(y,1),1)];

% linear tranformation
b = pinv(y)*x;
xfm = b;

% plot results
if ~nargout,
   figure;
   subplot(2,1,1);hold on;
   plot(x(:,1),x(:,2),'bo');
   plot(y(:,1),y(:,2),'rx');
   title('original data');
   legend('x','y');

   % interpolate
   yi=y*b;
   subplot(2,1,2);hold on;
   plot(x(:,1),x(:,2),'bo');
   plot(yi(:,1),yi(:,2),'rx');
   title('original data');
   legend('x','y interpolated');
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function spy( x, mode )
global cvx___

%   SPY is not a computational function, so it can be applied to any
%   CVX variables.

switch nargin,
    case 0,
        cvx_throw( 'Not enough arguments.' );
    case 1,
        mode = '';
    case 2,
        if ~ischar( mode ) || size( mode, 1 ) > 1,
            cvx_throw( 'Second argument must be a string.' );
        end
end

b = x.basis_;
s = x.size_;

switch mode,
    case { '2-d', '2-D', '2d', '2D' },
        b = sum( b ~= 0, 1 );
        b = reshape( b, s );
    case { '', '3-d', '3-D', '3d', '3D' },
        p = length( cvx___.classes ) + 1;
        if size( b, 1 ) < p, b( p, end ) = 0; end
    otherwise,
        cvx_throw( 'Unknown spy mode: %s', mode );
end

spy( b );

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

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


text:
function x = emailFeatures(word_indices)
%EMAILFEATURES takes in a word_indices vector and produces a feature vector
%from the word indices
%   x = EMAILFEATURES(word_indices) takes in a word_indices vector and
%   produces a feature vector from the word indices.

% Total number of words in the dictionary
n = 1899;

% You need to return the following variables correctly.
x = zeros(n, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return a feature vector for the
%               given email (word_indices). To help make it easier to
%               process the emails, we have have already pre-processed each
%               email and converted each word in the email into an index in
%               a fixed dictionary (of 1899 words). The variable
%               word_indices contains the list of indices of the words
%               which occur in one email.
%
%               Concretely, if an email has the text:
%
%                  The quick brown fox jumped over the lazy dog.
%
%               Then, the word_indices vector for this text might look
%               like:
%
%                   60  100   33   44   10     53  60  58   5
%
%               where, we have mapped each word onto a number, for example:
%
%                   the   -- 60
%                   quick -- 100
%                   ...
%
%              (note: the above numbers are just an example and are not the
%               actual mappings).
%
%              Your task is take one such word_indices vector and construct
%              a binary feature vector that indicates whether a particular
%              word occurs in the email. That is, x(i) = 1 when word i
%              is present in the email. Concretely, if the word 'the' (say,
%              index 60) appears in the email, then x(60) = 1. The feature
%              vector should look like:
%
%              x = [ 0 0 0 0 1 0 0 0 ... 0 0 0 0 1 ... 0 0 0 1 0 ..];
%
%

for i=1:length(word_indices)
    x(word_indices(i)) = 1;
end
% =========================================================================

end

author: merwan

repo: ml-class

save_path: github-repos/MATLAB/merwan-ml-class


text:
% op_alignISIS.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% [out,fs,phs]=op_alignISIS(in,tmax,initPars);
%
% DESCRIPTION:
% Apply spectral registration to align ISIS subspectra prior to
% subtraction.  This is intended to be used prior to averaging, so that the
% alignment can be performed independently for each average.
%
% INPUTS:
% in        = Input data structure.
% tmax      = Maximum time (s) in time domain to use for alignment.
% initPars	= (Optional) Initial fit parameters [freq(Hz), phase(degrees)]. Default=[0,0];
%
% OUTPUTS:
% out       = Output following alignment of ISIS subspectra.
% fs        = Vector of frequency shifts (in Hz) used for alignment.
% phs       = Vector of phase shifts (in degrees) used for alignment.

function [out,fs,phs]=op_alignISIS(in,tmax,initPars)

if ~in.flags.addedrcvrs
    error('ERROR:  I think it only makes sense to do this after you have combined the channels using op_addrcvrs.  ABORTING!!');
end

if in.dims.subSpecs==0
    error('ERROR:  Must have multiple subspectra.  ABORTING!!');
end

if nargin<3
    parsGuess=[0,0];
else
    parsGuess=initPars;
end

if in.dims.averages
    fs=zeros(in.sz(in.dims.averages),1);
    phs=zeros(in.sz(in.dims.averages),1);
else
    fs=0;
    phs=0;
end
fids=zeros(in.sz(in.dims.t),1);

disp('aligning all averages to the Average ISIS subtracted spectrum');
if in.dims.averages
    base0=op_median(op_combinesubspecs(in,'diff'));
else
    base0=op_combinesubspecs(in,'diff');
end
base=[real(base0.fids( base0.t>=0 & base0.t<tmax ));imag(base0.fids( base0.t>=0 & base0.t<tmax ))];
begin=1;
if in.dims.averages
    for n=begin:in.sz(in.dims.averages)
        %disp(['fitting subspec number ' num2str(m) ' and average number ' num2str(n)]);
        parsFit=nlinfit(squeeze(in.fids(in.t>=0 & in.t<tmax,n,:)),base,@op_freqPhaseShiftComplexNest,parsGuess);
        A=op_freqPhaseShiftNest(parsFit,in.fids(:,n,:));
        size(A);
        size(fids);
        fids(:,n,1)=A(:,1);
        fids(:,n,2)=A(:,2);
        fs(n)=parsFit(1);
        phs(n)=parsFit(2);
        %plot(in.ppm,fftshift(ifft(fids(:,1,m))),in.ppm,fftshift(ifft(fids(:,n,m))));
    end
else
    n=1;
    %disp(['fitting subspec number ' num2str(m) ' and average number ' num2str(n)]);
    parsFit=nlinfit(squeeze(in.fids(in.t>=0 & in.t<tmax,:)),base,@op_freqPhaseShiftComplexNest,parsGuess);
    A=op_freqPhaseShiftNest(parsFit,in.fids(:,:));
    size(A);
    size(fids);
    fids(:,1)=A(:,1);
    fids(:,2)=A(:,2);
    fs(n)=parsFit(1);
    phs(n)=parsFit(2);
    %plot(in.ppm,fftshift(ifft(fids(:,1,m))),in.ppm,fftshift(ifft(fids(:,n,m))));
end

%re-calculate Specs using fft
specs=fftshift(ifft(fids,[],in.dims.t),in.dims.t);

%FILLING IN DATA STRUCTURE
out=in;
out.fids=fids;
out.specs=specs;

%FILLING IN THE FLAGS
out.flags=in.flags;
out.flags.writtentostruct=1;
out.flags.freqcorrected=1;

    function y=op_freqPhaseShiftComplexNest(pars,input)
        f=pars(1);     %Frequency Shift [Hz]
        p=pars(2);     %Phase Shift [deg]

        dwelltime=in.dwelltime;
        t=0:dwelltime:(length(input)-1)*dwelltime;
        fid=input(:,2);

        shifted=addphase(fid.*exp(-1i*t'*f*2*pi),p);
        subtracted=input(:,1)+shifted;
        subtracted=subtracted/2;

        y=[real(subtracted);imag(subtracted)];
        %y=real(fid.*exp(-1i*t'*f*2*pi));

    end

    function y=op_freqPhaseShiftNest(pars,input)
        f=pars(1);     %Frequency Shift [Hz]
        p=pars(2);     %Phase Shift [deg]

        dwelltime=in.dwelltime;
        t=0:dwelltime:(length(input)-1)*dwelltime;
        fid=input(:,2);

        shifted=addphase(fid.*exp(-1i*t'*f*2*pi),p);
        y(:,1)=input(:,1);
        y(:,2)=shifted;
        %y=real(fid.*exp(-1i*t'*f*2*pi));

    end

end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function su_hdr = GE_readHeaderSuite(fid, byte_align)
%
% su_hdr = GE_readHeaderSuite(fid, byte_align)
%
% Loads the suite header from the file with filed id fid
% and returns it as a structure.
% if byte_align = 1 then 32-bit alignment (SGI, LX2 format)
% if byte_align = 0 then 16-bit alignment (Sun, 5.X format)
%
%
% Souheil J. Inati
% Dartmouth College
% May 2000
% souheil.inati@dartmouth.edu
%

% define the structure and read
su_hdr = struct('su_id', fread(fid,4,'char'));                   % Suite ID
su_hdr = setfield(su_hdr, 'su_uniq', fread(fid,1,'int16'));      % The Make-Unique Flag
su_hdr = setfield(su_hdr, 'su_diskid', fread(fid,1,'char'));     % Disk ID
su_hdr = setfield(su_hdr, 'prodid', fread(fid,13,'char'));       % Product ID
su_hdr = setfield(su_hdr, 'su_verscre', fread(fid,2,'char'));    % Genesis Version
su_hdr = setfield(su_hdr, 'su_verscur', fread(fid,2,'char'));    % Genesis Version
su_hdr = setfield(su_hdr, 'su_checksum', fread(fid,1,'uint32')); % Suite  Record Checksum
su_hdr = setfield(su_hdr, 'su_padding', fread(fid,85,'char'));   % Spare Space
fseek(fid,1,0); % 16-bit alignment
if byte_align, fseek(fid,2,0); end % 32-bit alignment

return


text:
function newquiver(X,Y,U,V,varargin)
% NEWQUIVER - Plot velocity vectors as arrows
%
% Syntax:
%   NEWQUIVER(X,Y,U,V,options)
%
% In:
%   X       - Grid of x coordinates
%   Y       - Grid of y coordinates
%   U       - Vector direction (x component)
%   V       - Vector direction (y component)
%   options - name-value pairs (see below)
%
% Out:
%   (none)
%
% Description:
%   NEWQUIVER(X,Y,U,V) plots velocity vectors as arrows with components
%   (u,v) at the points (x,y).  The matrices X,Y,U,V must all be the same
%   size and contain corresponding position and velocity components
%   (X and Y can also be vectors to specify a uniform grid).
%   NEWQUIVER automatically scales the arrows to fit within the grid.
%
%   Additional option-value pairs can be given.
%   Accepted options/values are:
%    * scale     : A two-element vector for scaling the shape
%    * color     : Fill color (rgbkcy...) (default 'k': black)
%    * EdgeColor : Edge color of shape  (default 'k': black)
%    * LineWidth : Line width (default 0.5)
%    * X         : Arrow shape as a two-column vector (default: arrow)
%
% See also:
%    quiver
%
% Copyright:
%   2012-2018 - Simo Särkkä and Arno Solin
%
% License:
%   This software is provided under the MIT License. See the accompanying
%   LICENSE file for details.

%% Define arrow shape

  % Define arrow shape
  arrow.X = [8.7185878,  4.0337352 % Upper right corner
            -2.2072895,  0         % Arrow tip
             8.7185878, -4.0337352 % Lower right corner
             6.9831476, -1.6157441 % Inner lower corner
             6.9831476, -0.8       % Inner lower shaft start
             28,        -0.8       % Lower shaft end
             28,         0.8       % Upper shaft end
             6.9831476,  0.8       % Inner upper shaft start
             6.9831476,  1.6157441 % Inner upper corner
             8.7185878,  4.0337352]; % Upper right corner

   % Flip arrow left to right
   arrow.X(:,1) = -arrow.X(:,1);

   % Normalize length to one and center
   arrow.X(:,1) = arrow.X(:,1)-min(arrow.X(:,1));
   arrow.X = arrow.X/max(arrow.X(:,1));
   arrow.X(:,1) = arrow.X(:,1)-.5;

%% Set default options

  % Default options structure
  defaultopt.scale       = .5*[min(X(X(:)>min(X(:))))-min(X(:)) ...
                               min(Y(Y(:)>min(Y(:))))-min(Y(:))];
  defaultopt.color       = 'k';
  defaultopt.EdgeColor   = 'k';
  defaultopt.LineWidth   = .5;
  defaultopt.X           = arrow.X;

%% Get options

  % If there are options submitted
  if length(varargin) > 1

    % Options
    options = cell2struct(varargin(2:2:end),varargin(1:2:end),2);

    % Merge the default and submitted options
    fnames = fieldnames(defaultopt);
    for i=1:length(fnames)
      %val = getfield(options,fnames{i});
      if ~isfield(options,fnames{i})
        options=setfield(options,fnames{i},getfield(defaultopt,fnames{i}));
      end
    end
  else
    options = defaultopt;
  end

  % Make sure there is a scaling factor in both directions
  if numel(options.scale)~=2, options.scale = [1 1]*options.scale(1); end

%% Get arrow directions

  % Calculate arrow direction
  TH = atan2(V,U);

  hold on

  for i=1:numel(X)

    % Roatate
    Z = [cos(TH(i)) -sin(TH(i)); sin(TH(i)) cos(TH(i))]*options.X';

    % Scale
    Z = bsxfun(@times,Z,options.scale(:));

    % Move
    Z = bsxfun(@plus,Z,[X(i); Y(i)]);

    % Plot
    patch(Z(1,:),Z(2,:),1,'FaceColor',options.color, ...
        'EdgeColor',options.EdgeColor, ...
        'LineWidth',options.LineWidth,'LineStyle','-')

  end

author: AaltoML

repo: SDE

save_path: github-repos/MATLAB/AaltoML-SDE


text:
function test_bug2511

% WALLTIME 00:20:00
% MEM 8gb
% DEPENDENCY ft_sourceplot ft_read_headshape

t1 = ft_read_mri(dccnpath('/home/common/matlab/spm8/canonical/single_subj_T1.nii'));
load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug2511.mat'));

% loading the data takes quite some time, as it is 4.7GB on disk, which is
% a bit silly, because only 1 variable is used, and the whole
% XXX.cfg.previous...... is not needed. JM reduced the file size on
% 20140629

cfg                  = [];
cfg.parameter        = 'avg.pow';
cfg.interpmethod    = 'sphere_avg';
sourceDiffIntAVG     = ft_sourceinterpolate(cfg, AverageSource ,t1);

cfgS                     = [];
cfgS.funparameter        = 'avg.pow';
cfgS.funcolormap         = 'jet';
cfgS.opacitylim          = [0 1];
cfgS.opacitymap          = 'rampup';
cfgS.method              = 'surface';
cfgS.projmethod          = 'nearest';
cfgS.interactive         = 'yes';
% cfgS.surffile            = 'surface_l4_both.mat'; % this was causing the problem, the file was removed from FT
cfgS.surffile            = 'surface_white_both.mat'; % this is the new surface that should be used instead
cfgS.surfdownsample      = 10;
cfgS.funcolorlim         = 'maxabs';
cfgS.interactive         = 'yes';
ft_sourceplot(cfgS, sourceDiffIntAVG);

%% error occurs with
% addpath(dccnpath('/home/common/matlab/fieldtrip/'));
%% but not with
% addpath(dccnpath('/home/common/matlab/fieldtrip-20131231/'));

%% Error description
%
% the input is volume data with dimensions [91 109 91]
% scaling anatomy to [0 1]
% no masking parameter
% The source data is defined on a 3D grid, interpolation to a surface mesh will be performed
% Warning: could not determine filetype of surface_l4_both.mat
% > In fileio/private/ft_warning at 158
%   In ft_filetype at 1189
%   In ft_read_headshape at 110
%   In ft_sourceplot at 890
% Error using ft_read_headshape (line 888)
% unknown fileformat "unknown" for head shape information
%
% Error in ft_sourceplot (line 890)
%     surf  = ft_read_headshape(cfg.surffile);


text:
function c=counts(x,v)
% COUNTS counts the number of elements of x that fall within v
% COUNTS(x,v) returns the number of elements such that
% v(1)<=x, ....v(n-1)<x<=v(n)
c=[];
nv=[-Inf; v];
for i=2:rows(nv);
   t=sum(x<=nv(i) & x>nv(i-1));
   c(i-1,1)=t;
end

author: Sable

repo: mcbench-benchmarks


text:
function cmap = polarmap(varargin)
%POLARMAP Polarized color map
%	POLARMAP applies a "polarized" blue-white-red colormap to current figure,
%	and adjusts the color axis limits to be centered to zero.
%
%	POLARMAP(M) fixes the number of colors to M (default is 64).
%
%	POLARMAP(MAP) applies linear shading to white to the center of colormap
%	MAP which can be any of existing colormaps (an Mx3 matrix of RGB).
%
%	POLARMAP(MAP,C) uses exponent C to modify the shading contrast. Default
%	is C = 1 for linear shading. Use C = 2 to strengthen the shading, or
%	C = 0.5 to attenuate it.
%
%	C=POLARMAP(...) returns an M-by-3 matrix containing the colormap, that
%	can be used with COLORMAP function like other colormaps.
%
%	Examples:
%		pcolor(peaks), shading interp
%		polarmap, colorbar
%
%	then try the following
%		polarmap(jet,0.5)
%
%	Note the polar shading has no real interest with colormaps that include
%	white color as one of the extremes (like GRAY, BONE, HOT, ...).
%
%	See also JET, HSV, COPPER, SPRING, SUMMER, WINTER, COOL, COLORMAP, RGBPLOT.
%
%	Author: Francois Beauducel, IPGP
%	Created: 2011-10-26
%	Updated: 2012-06-12

%	Copyright (c) 2012, François Beauducel, covered by BSD License.
%	All rights reserved.
%
%	Redistribution and use in source and binary forms, with or without
%	modification, are permitted provided that the following conditions are
%	met:
%
%	   * Redistributions of source code must retain the above copyright
%	     notice, this list of conditions and the following disclaimer.
%	   * Redistributions in binary form must reproduce the above copyright
%	     notice, this list of conditions and the following disclaimer in
%	     the documentation and/or other materials provided with the distribution
%
%	THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%	AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%	IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%	ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%	LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%	CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%	SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%	INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%	CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%	ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%	POSSIBILITY OF SUCH DAMAGE.

% default parameters
m = 64;	% number of colors
c = 1;	% exponent of shading factor (1 = linear)

if nargin > 0
	if ~isnumeric(varargin{1}) | (size(varargin{1},2) ~= 3 & ~isscalar(varargin{1}))
		error('First argument must be numeric: scalar M or Mx3 color matrix');
	end
	if isscalar(varargin{1})
		m = varargin{1};
	end
end
if nargin > 0 & size(varargin{1},2) == 3
		map = varargin{1};
		m = size(map,1);
else
	map = bluered(m);
end

if nargin > 1 & isscalar(varargin{2})
	c = varargin{2};
end

% linear shading from min/max (colormap value) to center (white)
r = repmat(abs(linspace(1,-1,m)).^c,[3,1])';
map = map.*r + 1 - r;

if nargout > 0
	cmap = map;
else
	colormap(map)
	caxis([-1,1]*max(abs(caxis)))
	% Note: this fixes color axis to manual mode...
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function map = bluered(m)

if mod(m,2)
	z = [0,0,0];
	m2 = floor(m/2);
else
	z = zeros([0,3]);
	m2 = m/2;
end
map = [repmat([0,0,1],[m2,1]);z;repmat([1,0,0],[m2,1])];


text:
function incilta() % autogenerated function wrapper
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals
    report_this_filefun(mfilename('fullpath'));
    %
    %
    % make another dialog interface and call cin_lta
    %
    %initial values
    iwl3 = 1;
    nustep = 10;

    figure(mess);
    clf
    set(gcf,'Name','Movie Parameter Input ');
    set(gca,'visible','off');
    set(gcf,'pos',[ ZG.welcome_pos 400 250]);

    % creates a dialog box to input some parameters
    %
    freq_field=uicontrol('Style','edit',...
        'Position',[.80 .60 .12 .10],...
        'Units','normalized','String',num2str(iwl3),...
        'callback',@callbackfun_001);

    inp2_field=uicontrol('Style','edit',...
        'Position',[.80 .40 .12 .10],...
        'Units','normalized','String',num2str(nustep),...
        'callback',@callbackfun_002);

    close_button=uicontrol('Style','Pushbutton',...
        'Position', [.60 .05 .15 .15 ],...
        'Units','normalized','callback',@callbackfun_003,'String','Cancel');

    go_button=uicontrol('Style','Pushbutton',...
        'Position',[.25 .05 .15 .15 ],...
        'Units','normalized',...
        'callback',@callbackfun_004,...
        'String','Go');

    txt1 = text(...
        'Position',[0. 0.65 0 ],...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Please input window length in years (iwl):');

    txt2 = text(...
        'Position',[0. 0.40 0 ],...
        'FontSize',ZmapGlobal.Data.fontsz.m ,...
        'FontWeight','bold',...
        'String','Please input number of movie-frames:');

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

    function callbackfun_001(mysrc,myevt)
        % automatically created callback function from text
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        iwl3=str2double(freq_field.String);
        freq_field.String=num2str(iwl3);
    end

    function callbackfun_002(mysrc,myevt)
        % automatically created callback function from text
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        nustep=str2double(inp2_field.String);
        inp2_field.String=num2str(nustep);
    end

    function callbackfun_003(mysrc,myevt)
        % automatically created callback function from text
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        zmapmenu;
    end

    function callbackfun_004(mysrc,myevt)
        % automatically created callback function from text
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        zmap_message_center();
        runcilta;
    end

end


text:
function [center,rad,v1n,v2nb] = circlefit3d(p1,p2,p3)
% circlefit3d: Compute center and radii of circles in 3d which are defined by three points on the circumference
% usage: [center,rad,v1,v2] = circlefit3d(p1,p2,p3)
%
% arguments: (input)
%  p1, p2, p3 - vectors of points (rowwise, size(p1) = [n 3])
%               describing the three corresponding points on the same circle.
%               p1, p2 and p3 must have the same length n.
%
% arguments: (output)
%  center - (nx3) matrix of center points for each triple of points in p1,  p2, p3
%
%  rad    - (nx1) vector of circle radii.
%           if there have been errors, radii is a negative scalar ( = error code)
%
%  v1, v2 - (nx3) perpendicular vectors inside circle plane
%
% Example usage:
%
%  (1)
%      p1 = rand(10,3);
%      p2 = rand(10,3);
%      p3 = rand(10,3);
%      [center, rad] = circlefit3d(p1,p2,p3);
%      % verification, result should be all (nearly) zero
%      result(:,1)=sqrt(sum((p1-center).^2,2))-rad;
%      result(:,2)=sqrt(sum((p2-center).^2,2))-rad;
%      result(:,3)=sqrt(sum((p3-center).^2,2))-rad;
%      if sum(sum(abs(result))) < 1e-12,
%       disp('All circles have been found correctly.');
%      else,
%       disp('There had been errors.');
%      end
%
%
% (2)
%       p1=rand(4,3);p2=rand(4,3);p3=rand(4,3);
%       [center,rad,v1,v2] = circlefit3d(p1,p2,p3);
%       plot3(p1(:,1),p1(:,2),p1(:,3),'bo');hold on;plot3(p2(:,1),p2(:,2),p2(:,3),'bo');plot3(p3(:,1),p3(:,2),p3(:,3),'bo');
%       for i=1:361,
%           a = i/180*pi;
%           x = center(:,1)+sin(a)*rad.*v1(:,1)+cos(a)*rad.*v2(:,1);
%           y = center(:,2)+sin(a)*rad.*v1(:,2)+cos(a)*rad.*v2(:,2);
%           z = center(:,3)+sin(a)*rad.*v1(:,3)+cos(a)*rad.*v2(:,3);
%           plot3(x,y,z,'r.');
%       end
%       axis equal;grid on;rotate3d on;
%
%
% Author: Johannes Korsawe
% E-mail: johannes.korsawe@volkswagen.de
% Release: 1.0
% Release date: 26/01/2012

% Default values
center = [];rad = 0;v1n=[];v2nb=[];

% check inputs
% check number of inputs
if nargin~=3,
    fprintf('??? Error using ==> cirlefit3d\nThree input matrices are needed.\n');rad = -1;return;
end
% check size of inputs
if size(p1,2)~=3 || size(p2,2)~=3 || size(p3,2)~=3,
    fprintf('??? Error using ==> cirlefit3d\nAll input matrices must have three columns.\n');rad = -2;return;
end
n = size(p1,1);
if size(p2,1)~=n || size(p3,1)~=n,
    fprintf('??? Error using ==> cirlefit3d\nAll input matrices must have the same number or rows.\n');rad = -3;return;
end
% more checks are to follow inside calculation

% Start calculation
% v1, v2 describe the vectors from p1 to p2 and p3, resp.
v1 = p2 - p1;v2 = p3 - p1;
% l1, l2 describe the lengths of those vectors
l1 = sqrt((v1(:,1).*v1(:,1)+v1(:,2).*v1(:,2)+v1(:,3).*v1(:,3)));
l2 = sqrt((v2(:,1).*v2(:,1)+v2(:,2).*v2(:,2)+v2(:,3).*v2(:,3)));
if find(l1==0) | find(l2==0), %#ok<OR2>
    fprintf('??? Error using ==> cirlefit3d\nCorresponding input points must not be identical.\n');rad = -4;return;
end
% v1n, v2n describe the normalized vectors v1 and v2
v1n = v1;for i=1:3, v1n(:,i) = v1n(:,i)./l1;end
v2n = v2;for i=1:3, v2n(:,i) = v2n(:,i)./l2;end
% nv describes the normal vector on the plane of the circle
nv = [v1n(:,2).*v2n(:,3) - v1n(:,3).*v2n(:,2) , v1n(:,3).*v2n(:,1) - v1n(:,1).*v2n(:,3) , v1n(:,1).*v2n(:,2) - v1n(:,2).*v2n(:,1)];
if find(sum(abs(nv),2)<1e-5),
    fprintf('??? Warning using ==> cirlefit3d\nSome corresponding input points are nearly collinear.\n');
end
% v2nb: orthogonalization of v2n against v1n
dotp = v2n(:,1).*v1n(:,1) + v2n(:,2).*v1n(:,2) + v2n(:,3).*v1n(:,3);
v2nb = v2n;for i=1:3,v2nb(:,i) = v2nb(:,i) - dotp.*v1n(:,i);end
% normalize v2nb
l2nb = sqrt((v2nb(:,1).*v2nb(:,1)+v2nb(:,2).*v2nb(:,2)+v2nb(:,3).*v2nb(:,3)));
for i=1:3, v2nb(:,i) = v2nb(:,i)./l2nb;end

% remark: the circle plane will now be discretized as follows
%
% origin: p1                    normal vector on plane: nv
% first coordinate vector: v1n  second coordinate vector: v2nb

% calculate 2d coordinates of points in each plane
% p1_2d = zeros(n,2); % set per construction
% p2_2d = zeros(n,2);p2_2d(:,1) = l1; % set per construction
p3_2d = zeros(n,2); % has to be calculated
for i = 1:3,
    p3_2d(:,1) = p3_2d(:,1) + v2(:,i).*v1n(:,i);
    p3_2d(:,2) = p3_2d(:,2) + v2(:,i).*v2nb(:,i);
end

% calculate the fitting circle
% due to the special construction of the 2d system this boils down to solving
% q1 = [0,0], q2 = [a,0], q3 = [b,c] (points on 2d circle)
% crossing perpendicular bisectors, s and t running indices:
% solve [a/2,s] = [b/2 + c*t, c/2 - b*t]
% solution t = (a-b)/(2*c)

a = l1;b = p3_2d(:,1);c = p3_2d(:,2);
t = 0.5*(a-b)./c;
scale1 = b/2 + c.*t;scale2 = c/2 - b.*t;

% centers
center = zeros(n,3);
for i=1:3,
    center(:,i) = p1(:,i) + scale1.*v1n(:,i) + scale2.*v2nb(:,i);
end

% radii
rad = sqrt((center(:,1)-p1(:,1)).^2+(center(:,2)-p1(:,2)).^2+(center(:,3)-p1(:,3)).^2);


text:
function [A,Phi] = inverse_cht(Am,Nphi)
%INVERSE_CHT inverse circular harmonics transform (ICHT)
%
%   Usage: [A,Phi] = inverse_cht(Am,[Nphi])
%
%   Input parameters:
%       Am      - circular harmonics coefficients [N x (2*M+1)]
%       Nphi    - number of equi-angular distributed angles, for which the ICHT
%                 is computed, optional, default: 2*M+1
%
%   Output parameters:
%       A       - inverse circular harmonics transform [N x Nphi]
%       Phi     - corresponding angle of the ICHT [1 x Nphi]
%
%   See also: pwd_imp_circexp

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 1;
nargmax = 2;
narginchk(nargmin,nargmax);
isargmatrix(Am);
if nargin == nargmin
    Nphi = size(Am, 2);
else
    isargpositivescalar(Nphi);
end

%% ===== Computation ==================================================
M = (size(Am,2)-1)/2;
N = size(Am,1);

% Implementation of
%           ___
%           \
% A(phi) =  /__    A  e^(+i*m*n*2*pi/Nphi)
%         m=-M..M   m

% Spatial IFFT
A = zeros(N, Nphi);
% this handles cases where Nphi < M
for l=1:N
    A(l,:) = sum(buffer(Am(l,:),Nphi),2);
end
A = circshift(A,[0,-M]);  % m = 0, ..., M, ..., -M, ..., -1
A = ifft(A,[],2) * Nphi;  % IFFT includes factor 1/Nphi

% Axis corresponding to ICHT
if nargout>1
    Phi = 0:2*pi / Nphi:2*pi*(1-1/Nphi);
end

author: sfstoolbox

repo: sfs-matlab


text:
function [ftps] = mph2ftps(mph)
% Convert speed from miles per hour to feet per second
% Chad A. Greene 2012
ftps = mph*1.466666666667;


text:
function b = bodewig_rhs ( )

%*****************************************************************************80
%
%% BODEWIG_RHS returns the BODEWIG right hand side.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real B(4,1), the right hand side vector.
%
  b(1:4,1) = [ 29.0, 18.0, 15.0, 4.0 ];

  return
end


text:
function success=testAddSecondExperimentToNewWingedFlyJab()

testDataDirName='/groups/branson/bransonlab/projects/JAABA/test_data/flies_wing';
%testDataDirName='/Users/taylora/jaaba/sampledata_copy/flies_wing';
nameOfExpToAdd1='pBDPGAL4U_TrpA_Rig2Plate17BowlA_20110929T143440';
nameOfExpToAdd2='GMR_71G01_AE_01_TrpA_Rig2Plate17BowlD_20110921T084818';

% Make a Macguffin, from which we'll open a new file in JLabelData
macguffin=Macguffin('flies_wing');
macguffin.setMainBehaviorName('may22ing');
macguffin.setScoreFileName('scores_may22ing.mat');
macguffin.setTrxFileName('wingtracking_results.mat');
macguffin.setMovieFileName('movie.ufmf');

% Create the JLabelData object
data=JLabelData('setstatusfn',@(str)(fprintf('%s\n',str)), ...
                'clearstatusfn',@()(nop()), ...
                'isInteractive',false);
data.perframeGenerate=true;  % tell JLabelData to generate any missing perframe files

% Create a new file in the JLabelData object
data.newJabFile(macguffin);

% Add a new exp dir
[success,msg]=data.AddExpDir(fullfile(testDataDirName,nameOfExpToAdd1));
if ~success ,
  error(msg);
end

% Add a second exp dir
[success,msg]=data.AddExpDir(fullfile(testDataDirName,nameOfExpToAdd2));
if ~success ,
  error(msg);
end

% Clean up
data.closeJabFile();
data=[];  %#ok
success=true;

end

author: kristinbranson

repo: JAABA


text:
% Construction of a bilinearly blended Coons surface
%

%  D.M. Spink
%  Copyright (c) 2000.

% Boundary curve 1
pnts = [ 0.0  3.0  4.5  6.5 8.0 10.0;
         0.0  0.0  0.0  0.0 0.0  0.0;
         2.0  2.0  7.0  4.0 7.0  9.0];
crv1 = nrbmak(pnts, [0 0 0 1/3 0.5 2/3 1 1 1]);

% Boundary curve 2
pnts= [ 0.0  3.0  5.0  8.0 10.0;
        10.0 10.0 10.0 10.0 10.0;
        3.0  5.0  8.0  6.0 10.0];
crv2 = nrbmak(pnts, [0 0 0 1/3 2/3 1 1 1]);

% Boundary curve 3
pnts= [ 0.0 0.0 0.0 0.0;
        0.0 3.0 8.0 10.0;
        2.0 0.0 5.0 3.0];
crv3 = nrbmak(pnts, [0 0 0 0.5 1 1 1]);

% Boundary curve 4
pnts= [ 10.0 10.0 10.0 10.0 10.0;
        0.0   3.0  5.0  8.0 10.0;
        9.0   7.0  7.0 10.0 10.0];
crv4 = nrbmak(pnts, [0 0 0 0.25 0.75 1 1 1]);

srf = nrbcoons(crv1, crv2, crv3, crv4);

% Draw the surface
nrbplot(srf,[20 20]);
title('Construction of a bilinearly blended Coons surface.');


text:
%% Generate best fit function
start_vy = [-0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8];
target_vy = [0];
start_vx = 0.0;%[-0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4];
target_vx = [0];
T = 0.4;
subfolder_name = 'library4';
if ~exist(fullfile('local', subfolder_name, 'transition_fitting'), 'dir')
    mkdir(fullfile('local', subfolder_name, 'transition_fitting'));
end

N_vx = length(start_vx);
N_vy = length(start_vy);
N = 16;
t_all = cell(N_vx, N_vy);
q_all = cell(N_vx, N_vy);
dq_all = cell(N_vx, N_vy);
px_all = cell(N_vx, N_vy);
py_all = cell(N_vx, N_vy);
vx_all = cell(N_vx, N_vy);
vy_all = cell(N_vx, N_vy);
counter = 1;
for i = 1:N_vx
    vx = start_vx(i);
    for j = 1:N_vy
        vy = start_vy(j);
        data_name = fullfile('local', subfolder_name, 'transition', ...
                sprintf('gait_X%0.1f_Y%.1f_TO_X%0.1f_Y%.1f.mat', vx, vy, target_vx, target_vy));
        param = load(data_name);

        t_all{i,j} = [param.gait(1).tspan]./T;%, param.gait(3).tspan];
        q_all{i,j} = [param.gait(1).states.x];%, param. gait(3).states.x];
        dq_all{i,j} = [param.gait(1).states.dx];%, param. gait(3).states.dx];
        px_all{i,j} = [param.gait(1).states.x(1,:)];%ones(size(t_all{i,j}))*vx;
        py_all{i,j} = [param.gait(1).states.x(2,:)];%ones(size(t_all{i,j}))*vy;
        vx_all{i,j} = [param.gait(1).states.dx(1,:)];%ones(size(t_all{i,j}))*vx;
        vy_all{i,j} = [param.gait(1).states.dx(2,:)];%ones(size(t_all{i,j}))*vy;

    end
end

inpt = [horzcat(t_all{:})
%     horzcat(px_all{:})
%     horzcat(py_all{:})
%     horzcat(vx_all{:})
    horzcat(vy_all{:})];
q_outpt = horzcat(q_all{:});
dq_outpt = horzcat(dq_all{:});
%%
for i = 1:N
    q_outpt_specific = q_outpt(i, :);
    dq_outpt_specific = dq_outpt(i, :);
    opt.neuralFitting(q_outpt_specific, inpt, 5, fullfile('local', subfolder_name, 'transition_fitting', sprintf('sagittal_library_q%i', i)));
    opt.neuralFitting(dq_outpt_specific, inpt, 5, fullfile('local', subfolder_name, 'transition_fitting', sprintf('sagittal_library_dq%i', i)));
end
addpath(fullfile('local', subfolder_name, 'transition_fitting'));
%% Plot features
% joint angle
joint_names = {    'BasePosX'
    'BasePosY'
    'BasePosZ'
    'BaseRotX'
    'BaseRotY'
    'BaseRotZ'
    'qHRight'
    'qARight'
    'qBRight'
    'fourBarARight'
    'fourBarBRight'
    'qHLeft'
    'qALeft'
    'qBLeft'
    'fourBarALeft'
    'fourBarBLeft'};
for n = 1:N
    f = figure(1000+n);
    f.Name = joint_names{n};
    set(f, 'WindowStyle', 'docked');

    ax = axes(f);
    hold(ax);

    scatter3(ax, inpt(1,:),inpt(2,:),q_outpt(n,:), 'r');

    [S, V] = meshgrid(0:0.01:1, -2:0.04:2);
    L = zeros(size(S));
    for i = 1:size(S, 1)
        for j = 1:size(S, 2)
            L(i, j) = feval(sprintf('sagittal_library_q%i', n), ([S(i, j); V(i, j)]));
            %L(i, j) = l(n);
        end
    end

    surface(ax, S, V, L);
end
%% joint velocity
for n = 1:N
    f = figure(2000+n);
    f.Name = joint_names{n};
    set(f, 'WindowStyle', 'docked');

    ax = axes(f);
    hold(ax);

    scatter3(ax, inpt(1,:),inpt(2,:),dq_outpt(n,:), 'r');

    [S, V] = meshgrid(0:0.01:1, -2:0.01:2);
    L = zeros(size(S));
    for i = 1:size(S, 1)
        for j = 1:size(S, 2)
            L(i, j) = feval(sprintf('sagittal_library_dq%i', n), ([S(i, j); V(i, j)]));
        end
    end

    surface(ax, S, V, L);
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function cl = orthviews(image_obj, varargin)
% Orthviews display (SPM) for CANlab image_vector (or fmri_data, statistic_image) object
%
% *Usage:
% ::
%
%    orthviews(image_obj, varargin)
%
% :Optional Inputs:
%   **'posneg':**
%        input generates orthviews using solid colors, separated for positive- and negative-valued voxels.
%
%   **'largest_region':**
%        to center the orthviews on the largest region in the image
%
%   **'overlay':**
%        followed by name of anatomical image to use as underlay
%
%   **'unique':**
%        plot groups of contiguous voxels in different, unique colors
%        - autodetect feature turns this on when there are <300 unique
%        values, or for 'atlas' objects. See 'continuous' below to force continuous.
%
%   **'parcels':**
%        plot voxels with each unique value in a different color
%        useful for images in which values indicate parcel or network
%        membership
%        Default is 'unique' for <300 values, 'continuous' for more.  Use
%        'unique' to force unique-valued colors.
%
%   **'continuous':**
%        Plot voxels color-mapped according to their values.
%        Default is 'unique' for <300 values, 'continuous' for more.  Use
%        'unique' to force unique-valued colors.
%
%   **'input_handle':**
%        followd by integer for which orthviews display to add to.
%        spm_orthviews can display multiple linked orthviews sections. If
%        you have initialized a display with multiple images, Enter 1 to
%        add to the first one, 2 for the second, and so forth.
%
%   **'add':**
%        Add to the first orthviews display (e.g., 'input_handle', 1)
%
%   **'color':**
%       Followed by cell with 3-element rgb vector,
%       e.g., 'color', {[1 0 1]}
%       This is superseded by 'unique' and 'posneg' options.
%
% ..
%    Copyright Tor Wager, 2011
% ..

% Aug 2017: Tor Wager edited to allow flexibility in colors; unique,
% parcel, continuous. and autodetect.

%overlay = which('SPM8_colin27T1_seg.img');
overlay = which('keuken_2014_enhanced_for_underlay.img');

doposneg = false;
doreg = false;
input_handle=[];
dounique = false;
doparcels = false;
force_continuous = false;
additional_inputs = {};

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}
            % reserved keywords
            case 'posneg',  doposneg = 1;

                % functional commands
            case 'overlay', overlay = varargin{i + 1}; varargin{i + 1} = [];

            case {'han', 'handle', 'input_handle'}, input_handle = varargin{i+1};

            case 'add', input_handle = 1;

            case {'largest_region', 'largest_cluster'}, doreg = 1;

            case 'unique'
                dounique = true;

            case 'continuous'
                force_continuous = true;

            case 'parcels'
                doparcels = true; dounique = true;

            case 'color'
                additional_inputs{end + 1} = varargin{i};

                if ~iscell(varargin{i + 1})
                    warning('Color input must be a cell array. Check inputs');
                    varargin{i + 1} = {varargin{i + 1}};
                end

                additional_inputs{end + 1} = varargin{i + 1};

            case 'trans'
                additional_inputs{end + 1} = 'trans';

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

if ~exist(overlay, 'file')
    overlay = which(overlay);
end

if isempty(overlay) || ~exist(overlay, 'file')
    disp('Cannot find overlay image:')
    disp(overlay)
    return
end

if isempty(input_handle)
    % re-initialize new orthviews

    n = size(image_obj.dat, 2);
    if n > 24
        disp('Warning! Only first 24 images can displayed with spm_orthviews.');
        n = 24;
    end

    spm_check_registration(repmat(overlay, n, 1));

    handle_indices = 1:n;

else

    % use existing
    handle_indices = input_handle;

end

% replace missing voxels if necessary
image_obj = replace_empty(image_obj);

% enforce double, same var types across all objects
image_obj.dat = double(image_obj.dat);

% autodetect parcels and turn on 'unique' for maps with few values
% or atlas objects.
% (default = do this)
if isa(image_obj, 'atlas') && ~force_continuous
    dounique = true; doparcels = true;
end

u = unique(image_obj.dat(:));

if ~doparcels && ~force_continuous && all(u == round(u))
    % auto-select parcel unique color plot if there are few unique values and values are integers
    nvalues = length(u);

    if nvalues > 2 && nvalues < 300, dounique = true; doparcels = true; end
end

for i = handle_indices

    wh_image = min(i, size(image_obj.dat,2));  %min() is to differentiate between which subplot to plot on and which image to plot.  used when this is called with orthviews_multiple_objs.  Yoni 11/14

    if doparcels

        cl{i} = region(get_wh_image(image_obj, wh_image), 'unique_mask_values');

    else
        cl{i} = region(get_wh_image(image_obj, wh_image), 'contiguous_regions');

    end

    % old:
    %cl{i} = iimg_indx2clusters(image_obj.dat(:, wh_image), image_obj.volInfo);

    if dounique

        %colors = custom_colors([1 1 0], [0 0 1], length(cl{i}));

        if ~isempty(cl{i})
            cluster_orthviews(cl{i}, 'unique', 'add', 'handle', i, additional_inputs{:});
        else
            fprintf('Image %3.0f empty\n', i);
        end

    elseif doposneg

        warning off
        cl{i} = cluster2region(cl{i});
        warning on % name field warning

        [clpos{i}, clneg{i}] = posneg_separate(cl{i}, 'Z');

        if ~isempty(clpos{i})
            cluster_orthviews(clpos{i}, {[1 1 0]}, 'add', 'handle', i,  'solid'); % additional_inputs{:}); % 'solid');
        else
            fprintf('Image %3.0f: No positive clusters\n', i);
        end

        if ~isempty(clneg{i})
            cluster_orthviews(clneg{i}, {[0 0 1]}, 'add', 'handle', i,  'solid'); % additional_inputs{:}); % 'solid');
        else
            fprintf('Image %3.0f: No negative clusters\n', i);
        end

    else
        if ~isempty(cl{i})
            cluster_orthviews(cl{i}, additional_inputs{:}, 'add', 'handle', i, additional_inputs{:});
        else
            fprintf('Image %3.0f empty\n', i);
        end
    end

end

% Set colormap
% If one image, we can use 'hot/cool' split, which shows zero point
% If multiple images, default to a single colormap, as the entire
% figure can use only one colormap. Note: colormaps can apply to specific
% axes...

if length(handle_indices) > 1
    %spm_orthviews_change_colormap([0 0 1], [1 1 0], [0 1 1], [.5 .5 .5], [1 .5 0]);
    %spm_orthviews_change_colormap([.5 0 1], [1 1 0]);
    spm_orthviews_change_colormap([.2 .2 .6], [1 1 0]);  % slate to yellow

    % won't work now -
    %     maxposcolor = [1 1 0]; % max pos, most extreme values
    %     minposcolor = [1 .4 .5]; % [.8 .3 0]; % min pos
    %     maxnegcolor = [0 .8 .8]; % [.3 .6 .9]; % max neg
    %     minnegcolor = [0 0 1]; % min neg, most extreme values
    %
    %     spm_orthviews_change_colormap(minnegcolor, maxnegcolor, minposcolor, maxposcolor);

else % one image.
    spm_orthviews_hotcool_colormap(image_obj.dat(:), 0);
end

% center image
spm_orthviews('Reposition', [0 0 0]);

% if requested, try to center on largest region
if doreg
    [~,wh] = max(cat(1,cl{1}.numVox));
    spm_orthviews('Reposition', cl{1}(wh).mm_center);
end
drawnow;

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [hmm,Gamma,vpath,Xi,LL,datat] = hmmdual(data,T,hmm,Gamma,Xi,residuals)
%
% Dual estimation of the HMM, first Gamma and then the HMM structure
%
% INPUTS:
%
% data          observations - a struct with X (time series) and C (classes; optional)
% T             Number of time points for each time series
% hmm           hmm structure with options specified in hmm.train
% Gamma         Initial state courses
% Xi            joint probability of past and future states conditioned on data
% residuals     in case we train on residuals, the value of those.
%
% OUTPUTS
% hmm           estimated HMMMAR model
% Gamma         estimated p(state | data)
% vpath            estimated Viterbi path
%
% Author: Diego Vidaurre, OHBA, University of Oxford (2019)
%
% edits to work for gradient computation (for Fisher kernel):
% Christine Ahrends, Aarhus University 2022

% to fix potential compatibility issues with previous versions
hmm = versCompatibilityFix(hmm);

if nargin<6, residuals = []; end
if nargin<5, Xi = []; end
if nargin<4, Gamma = []; end

if iscell(T)
    for i = 1:length(T)
        if size(T{i},1)==1, T{i} = T{i}'; end
    end
    if size(T,1)==1, T = T'; end
    T = cell2mat(T);
end
checkdatacell;
N = length(T);
p = hmm.train.lowrank; do_HMM_pca = (p > 0);

train = hmm.train;
checkdatacell;
data = data2struct(data,T,train);
% if train.standardise
%     disp('Option standardise should be zero in hmmdual.')
%     disp('Standardization should be done separately and using the entire data set.')
% end
% Standardise data and control for ackward trials
valid_dims = computeValidDimensions(data,train);
data = standardisedata(data,T,train.standardise,valid_dims);
% Filtering
if ~isempty(train.filter)
    data = filterdata(data,T,train.Fs,train.filter);
end
% Detrend data
if train.detrend
    data = detrenddata(data,T);
end
% Leakage correction
if train.leakagecorr ~= 0
    data = leakcorr(data,T,train.leakagecorr);
end
% Hilbert envelope
if train.onpower
    data = rawsignal2power(data,T);
end
% Leading Phase Eigenvectors
if train.leida
    data = leadingPhEigenvector(data,T);
end
% pre-embedded  PCA transform
if length(train.pca_spatial) > 1 || train.pca_spatial > 0
    if isfield(train,'As')
        data.X = bsxfun(@minus,data.X,mean(data.X));
        data.X = data.X * train.As;
    else
        [train.As,data.X] = highdim_pca(data.X,T,train.pca_spatial);
    end
end
% Embedding
if length(train.embeddedlags) > 1
    [data,T] = embeddata(data,T,train.embeddedlags);
end
% PCA transform
if length(train.pca) > 1 || train.pca > 0
    if isfield(train,'A')
        data.X = bsxfun(@minus,data.X,mean(data.X));
        data.X = data.X * train.A;
    else
        error('PCA cannot be recomputed within hmmdual, use parameter A instead')
        %[train.A,data.X] = highdim_pca(data.X,T,train.pca,0,0,0,train.varimax);
    end
    % Standardise principal components and control for ackward trials
    data = standardisedata(data,T,train.standardise_pc);
    train.ndim = size(train.A,2);
    train.S = ones(train.ndim);
    orders = formorders(train.order,train.orderoffset,train.timelag,train.exptimelag);
    train.Sind = formindexes(orders,train.S) == 1;
end
% Downsampling
if train.downsample > 0
    [data,T] = downsampledata(data,T,train.downsample,train.Fs);
end

if isempty(residuals) && ~do_HMM_pca
    if ~isfield(hmm.train,'Sind')
        orders = formorders(hmm.train.order,hmm.train.orderoffset,hmm.train.timelag,hmm.train.exptimelag);
        hmm.train.Sind = formindexes(orders,hmm.train.S) == 1;
    end
    residuals =  getresiduals(data.X,T,hmm.train.S,hmm.train.maxorder,hmm.train.order,...
        hmm.train.orderoffset,hmm.train.timelag,hmm.train.exptimelag,hmm.train.zeromean);
end

if isempty(Gamma)
    [Gamma,~,Xi,LL] = hsinference(data,T,hmm,residuals);
elseif isempty(Xi)
    Xi = approximateXi(Gamma,T,hmm);
end
setxx;

hmm = obsupdate(Gamma,hmm,residuals,XX,XXGXX);
hmm = hsupdate(Xi,Gamma,T,hmm);

if nargout > 1
    Gamma = hsinference(data,T,hmm,residuals);
end
if nargout > 2
    vpath = hmmdecode(data,T,hmm,1,residuals,0);
end
datat = data;

end

author: OHBA-analysis

repo: HMM-MAR


text:
function [energy,constrast,Entropy,Homogeneity,standard_dev,Ph,Slope] = getHaralicParams(structNum,planC)
%function getHaralicParams(structNum)
%
%This function returns the haralic parameters for structure structNum.
%
%APA,12/26/2006

if ~exist('planC')
    global planC
end
indexS = planC{end};

scanNum                             = getStructureAssociatedScan(structNum,planC);
[rasterSegments, planC, isError]    = getRasterSegments(structNum,planC);
[mask3M, uniqueSlices]              = rasterToMask(rasterSegments, scanNum, planC);
scanArray3M                         = getScanArray(planC{indexS.scan}(scanNum));
SUVvals3M                           = mask3M.*double(scanArray3M(:,:,uniqueSlices));
[minr, maxr, minc, maxc, mins, maxs]= compute_boundingbox(mask3M);
volToEval                           = SUVvals3M(minr:maxr,minc:maxc,mins:maxs);
volToEval(volToEval==0)             = NaN;
%volToEval                           = volToEval - min(volToEval(:));
volToEval                           = volToEval / max(volToEval(:));
%volToEval                           = sqrt(volToEval);
[f,Ph]                              = haralick3D(volToEval,16);
suv3M                               = scanArray3M(:,:,uniqueSlices);
maskScaled3D = suv3M(find(mask3M));
% maskScaled3D = maskScaled3D/
standard_dev                        = std(single(maskScaled3D));
%standard_dev                        = std(suv3M(find(mask3M)));
energy                              = f(1);
constrast                           = f(2);
Entropy                             = f(4);
Homogeneity                         = f(8);
%Calculate slope
indexS = planC{end};
init_th  = 10;
final_th = 80;
n_th     = 10;
Thresholds = linspace (init_th, final_th, n_th);
% Slope = calc_slope_grigsby(structNum,Thresholds,planC);
Slope = NaN;
disp(['Energy     : ',num2str(energy)])
disp(['Contrast   : ',num2str(constrast)])
disp(['Entropy    : ',num2str(Entropy)])
disp(['Homogenity : ',num2str(Homogeneity)])
disp(['Std. Dev.  : ',num2str(standard_dev)])
disp(['Slope      : ',num2str(Slope)])

author: cerr

repo: CERR


text:
function X = imNormalize( X, flag )
% Various ways to normalize a (multidimensional) image.
%
% X may have arbitrary dimension (ie an image or video, etc).  X is treated
% as a vector of pixel values.  Hence, the mean of X is the average pixel
% value, and likewise the standard deviation is the std of the pixels from
% the mean pixel.
%
% USAGE
%  X = imNormalize( X, flag )
%
% INPUTS
%  X       - n dimensional array to standardize
%  flag    - [1] determines normalization procedure. Sets X to:
%            1: have zero mean and unit variance
%            2: range in [0,1]
%            3: have zero mean
%            4: have zero mean and unit magnitude
%            5: zero mean/unit variance, throws out extreme values
%               and also normalizes to [0,1]
%
% OUTPUTS
%  X       - X after normalization.
%
% EXAMPLE
%  I = double(imread('cameraman.tif'));
%  N = imNormalize(I,1);
%  [mean(I(:)), std(I(:)), mean(N(:)), std(N(:))]
%
% See also FEVALARRAYS
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

if (isa(X,'uint8')); X = double(X); end
if (nargin<2 || isempty(flag)); flag=1; end
siz = size(X);

if( flag==1 || flag==3 || flag==4 )
  % set X to have zero mean
  X = X(:);  n = length(X);
  meanX = sum(X)/n;
  X = X - meanX;

  % set X to have unit std
  if( flag==1 || flag==4 )
    sumX2 = sum(X.^2);
    if( sumX2>0 )
      if( flag==4 )
        X = X / sqrt(sumX2);
      else
        X = X / sqrt(sumX2/n);
      end
    end
  end
  X = reshape(X,siz);

elseif(flag==2)
  % set X to range in [0,1]
  X = X - min(X(:));  X = X / max(X(:));

elseif( flag==5 )
  X = imNormalize( X, 1 );
  t=2;
  X( X<-t )= -t;
  X( X >t )=  t;
  X = X/2/t + .5;

else
  error('Unknown standardization procedure');
end

author: zouchuhang

repo: LayoutNet


text:
function [x,P]=ukf(fstate,x,P,hmeas,z,Q,R)
% UKF   Unscented Kalman Filter for nonlinear dynamic systems
% [x, P] = ukf(f,x,P,h,z,Q,R) returns state estimate, x and state covariance, P
% for nonlinear dynamic system (for simplicity, noises are assumed as additive):
%           x_k+1 = f(x_k) + w_k
%           z_k   = h(x_k) + v_k
% where w ~ N(0,Q) meaning w is gaussian noise with covariance Q
%       v ~ N(0,R) meaning v is gaussian noise with covariance R
% Inputs:   f: function handle for f(x)
%           x: "a priori" state estimate
%           P: "a priori" estimated state covariance
%           h: fanction handle for h(x)
%           z: current measurement
%           Q: process noise covariance
%           R: measurement noise covariance
% Output:   x: "a posteriori" state estimate
%           P: "a posteriori" state covariance
%
% Example:
%{
n=3;      %number of state
q=0.1;    %std of process
r=0.1;    %std of measurement
Q=q^2*eye(n); % covariance of process
R=r^2;        % covariance of measurement
f=@(x)[x(2);x(3);0.05*x(1)*(x(2)+x(3))];  % nonlinear state equations
h=@(x)x(1);                               % measurement equation
s=[0;0;1];                                % initial state
x=s+q*randn(3,1); %initial state          % initial state with noise
P = eye(n);                               % initial state covraiance
N=20;                                     % total dynamic steps
xV = zeros(n,N);          %estmate        % allocate memory
sV = zeros(n,N);          %actual
zV = zeros(1,N);
for k=1:N
  z = h(s) + r*randn;                     % measurments
  sV(:,k)= s;                             % save actual state
  zV(k)  = z;                             % save measurment
  [x, P] = ukf(f,x,P,h,z,Q,R);            % ekf
  xV(:,k) = x;                            % save estimate
  s = f(s) + q*randn(3,1);                % update process
end
for k=1:3                                 % plot results
  subplot(3,1,k)
  plot(1:N, sV(k,:), '-', 1:N, xV(k,:), '--')
end
%}
% Reference: Julier, SJ. and Uhlmann, J.K., Unscented Filtering and
% Nonlinear Estimation, Proceedings of the IEEE, Vol. 92, No. 3,
% pp.401-422, 2004.
%
% By Yi Cao at Cranfield University, 04/01/2008
%
L=numel(x);                                 %numer of states
m=numel(z);                                 %numer of measurements
alpha=1e-3;                                 %default, tunable
ki=0;                                       %default, tunable
beta=2;                                     %default, tunable
lambda=alpha^2*(L+ki)-L;                    %scaling factor
c=L+lambda;                                 %scaling factor
Wm=[lambda/c 0.5/c+zeros(1,2*L)];           %weights for means
Wc=Wm;
Wc(1)=Wc(1)+(1-alpha^2+beta);               %weights for covariance
c=sqrt(c);
X=sigmas(x,P,c);                            %sigma points around x
[x1,X1,P1,X2]=ut(fstate,X,Wm,Wc,L,Q);          %unscented transformation of process
% X1=sigmas(x1,P1,c);                         %sigma points around x1
% X2=X1-x1(:,ones(1,size(X1,2)));             %deviation of X1
[z1,Z1,P2,Z2]=ut(hmeas,X1,Wm,Wc,m,R);       %unscented transformation of measurments
P12=X2*diag(Wc)*Z2';                        %transformed cross-covariance
K=P12*inv(P2);
x=x1+K*(z-z1);                              %state update
P=P1-K*P12';                                %covariance update

function [y,Y,P,Y1]=ut(f,X,Wm,Wc,n,R)
%Unscented Transformation
%Input:
%        f: nonlinear map
%        X: sigma points
%       Wm: weights for mean
%       Wc: weights for covraiance
%        n: numer of outputs of f
%        R: additive covariance
%Output:
%        y: transformed mean
%        Y: transformed smapling points
%        P: transformed covariance
%       Y1: transformed deviations

L=size(X,2);
y=zeros(n,1);
Y=zeros(n,L);
for k=1:L
    Y(:,k)=f(X(:,k));
    y=y+Wm(k)*Y(:,k);
end
Y1=Y-y(:,ones(1,L));
P=Y1*diag(Wc)*Y1'+R;

function X=sigmas(x,P,c)
%Sigma points around reference point
%Inputs:
%       x: reference point
%       P: covariance
%       c: coefficient
%Output:
%       X: Sigma points

A = c*chol(P)';
Y = x(:,ones(1,numel(x)));
X = [x Y+A Y-A];


text:
function [nvar, names, types] = nex_info(filename)
% nex_info(filename) -- read and display .nex file info
%
% [nvar, names, types] = nex_info(filename)
%
% INPUT:
%   filename - if empty string, will use File Open dialog
% OUTPUT:
%   nvar - number of variables in the file
%   names - [nvar 64] array of variable names
%   types - [1 nvar] array of variable types
%           Interpretation of type values: 0-neuron, 1-event, 2-interval, 3-waveform,
%                        4-population vector, 5-continuous variable, 6 - marker

if(nargin ~= 1)
   disp('1 input arguments are required')
   return
end

if(isempty(filename))
   [fname, pathname] = uigetfile('*.nex', 'Select a Nex file');
	filename = strcat(pathname, fname);
end

fid = fopen(filename, 'r');
if(fid == -1)
	disp('cannot open file');
   return
end

disp(strcat('file = ', filename));
magic = fread(fid, 1, 'int32');
version = fread(fid, 1, 'int32');
comment = fread(fid, 256, 'char');
freq = fread(fid, 1, 'double');
tbeg = fread(fid, 1, 'int32');
tend = fread(fid, 1, 'int32');
nvar = fread(fid, 1, 'int32');
fseek(fid, 260, 'cof');
disp(strcat('version = ', num2str(version)));
disp(strcat('frequency = ', num2str(freq)));
disp(strcat('duration (sec) = ', num2str((tend - tbeg)/freq)));
disp(strcat('number of variables = ', num2str(nvar)));
names = zeros(1, 64);
types=zeros(nvar);
for i=1:nvar
	types(i) = fread(fid, 1, 'int32');
	var_version = fread(fid, 1, 'int32');
	names(i, :) = fread(fid, [1 64], 'char');
	dummy = fread(fid, 128+8, 'char');
end
names = char(names);
fclose(fid);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [fdx, fdy, fdz] = calcFdr(N, directionFlag)
%CALCFDR Function who's
%diagonals are the k-space representation of the image-domain differencing
%operator. See "vx" in the reference below.
%
%   N: Size of volume
%   directionFlag: 'forward' or 'backward', direction of the
%   differentiation.
%
%   Code refractored from Berkin Bilgic's scripts: "script_Laplacian_unwrap_Sharp_Fast_TV_gre3D.m"
%   and "script_Laplacian_unwrap_Sharp_Fast_TV_gre3D.m"
%   Original source: https://martinos.org/~berkin/software.html
%
%   Original reference:
%   Bilgic et al. (2014), Fast quantitative susceptibility mapping with
%   L1-regularization and automatic parameter selection. Magn. Reson. Med.,
%   72: 1444-1459. doi:10.1002/mrm.25029
%

    [k2,k1,k3] = meshgrid(0:N(2)-1, 0:N(1)-1, 0:N(3)-1);

    switch directionFlag
        case 'forward'
            fdx = 1 - exp(-2*pi*1i*k1/N(1));
            fdy = 1 - exp(-2*pi*1i*k2/N(2));
            fdz = 1 - exp(-2*pi*1i*k3/N(3));
        case 'backward'
            fdx = -1 + exp(2*pi*1i*k1/N(1));
            fdy = -1 + exp(2*pi*1i*k2/N(2));
            fdz = -1 + exp(2*pi*1i*k3/N(3));
    end

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [distAll, scoresAll, pointsAll] = assignGTinference(keypointsAll,annolist,pidxsAll,parts,thresh_gt,marg_scores,nTrials)

if (nargin < 5)
    thresh_gt = 0.5;
end

scoresAll = cell(length(pidxsAll),1);
pointsAll = cell(length(pidxsAll),1);
distAll = cell(length(pidxsAll),1);

for i = 1:length(pidxsAll)
    scoresAll{i} = cell(length(annolist),1);
    pointsAll{i} = cell(length(annolist),1);
    distAll{i} = cell(length(annolist),1);
end

for imgidx = 1:length(annolist)

%     figure(200); clf; imagesc(imread(annolist(imgidx).image.name));
%     hold on; axis equal;

    dist = cell(2,1);

    det_c = {'g','b'};
    for detidx = 1:nTrials

        dist{detidx} = inf(length(pidxsAll),length(annolist(imgidx).annorect));
        for i = 1:length(pidxsAll)
            pidx = pidxsAll(i);
            % part is a joint
            assert(parts(pidx+1).pos(1) == parts(pidx+1).pos(2));
            jidx = parts(pidx+1).pos(1);

            det = keypointsAll(imgidx).det{jidx+1};
            if (~isempty(det))
                [val,id] = max(det(:,4+detidx));
                pp = det(id,1:2);
            else
                val = -inf;
                pp = [0 0];
            end
            pointsAll{i}{imgidx} = [pointsAll{i}{imgidx}; pp];
            if (marg_scores)
                scoresAll{i}{imgidx} = [scoresAll{i}{imgidx}; val];
            else
                scoresAll{i}{imgidx} = [scoresAll{i}{imgidx}; det(id,3)];
            end

            for ridx = 1:length(annolist(imgidx).annorect)
                rect = annolist(imgidx).annorect(ridx);
                refDist = util_get_head_size(rect);
                if (isfield(rect, 'annopoints') && isfield(rect.annopoints, 'point'))
                    p = util_get_annopoint_by_id(rect.annopoints.point, jidx);
                    if (~isempty(p))
                        dist{detidx}(i,ridx) = norm([p.x p.y] - pp)/refDist;
%                         plot(p.x,p.y,'ro','MarkerFaceColor','y','MarkerEdgeColor','k','MarkerSize',10);
                    else
                        dist{detidx}(i,ridx) = nan;
                    end
                end
            end
%             plot(pp(1),pp(2),'ro','MarkerFaceColor',det_c{detidx},'MarkerEdgeColor','k','MarkerSize',10);
        end
    end

    for detidx = 1:2
        pckSum = zeros(1,size(dist{detidx},2));
        idxsNotNaN = ~isnan(dist{detidx});
        for i=1:size(dist{detidx},2)
            pckSum(i) = sum(dist{detidx}(idxsNotNaN(:,i),i) <= thresh_gt);
        end
        [val,idx] = max(pckSum);
        idx2 = setdiff(1:size(dist{detidx},2),idx);
        dist{detidx}(idxsNotNaN(:,idx2),idx2) = inf;
        for i = 1:length(pidxsAll)
            distAll{i}{imgidx} = [distAll{i}{imgidx}; dist{detidx}(i,:)];
        end
    end

end

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
dbclear all;

% compile matlab wrappers
disp('Building wrappers ...');
mex('readTrackletsMex.cpp','-I../cpp','-lboost_serialization');
disp('...done!');

repo: MSCKF


text:
function opts = mergeOptions(opts1, opts2)
% Merges two options structures with one having precedence over the other.
%
% function opts = mergeOptions(opts1, opts2)
%
% input: opts1 and opts2 are two structures.
% output: opts is a structure containing all fields of opts1 and opts2.
% Whenever a field is present in both opts1 and opts2, it is the value in
% opts2 that is kept.
%
% The typical usage is to have opts1 contain default options and opts2
% contain user-specified options that overwrite the defaults.
%
% See also: getGlobalDefaults

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

    if isempty(opts1)
        opts1 = struct();
    end
    if isempty(opts2)
        opts2 = struct();
    end

    opts = opts1;
    fields = fieldnames(opts2);
    for i = 1 : length(fields)
        opts.(fields{i}) = opts2.(fields{i});
    end

end

author: andrewssobral

repo: lrslibrary


text:
function [Nbegcol,Nlencol,Nrowndx,Nobjcnt,Nobjndx,cJacobian] = jacSparsity(interfacedata)

linear = setdiff(find(interfacedata.variabletype == 0),interfacedata.evalVariables);
oJacobian = zeros(length(linear),1);
variabletype = interfacedata.variabletype;
c = interfacedata.c;
F_struc = interfacedata.F_struc;
m = size(interfacedata.F_struc,1);

nonlinear = variabletype;
nonlinear(interfacedata.evalVariables) = 1;

for i = find(c(:)')%1:length(c)
%    if c(i)
        if nonlinear(i)%variabletype(i) | ismember(i,interfacedata.evalVariables)
            if ismember(i,interfacedata.evalVariables)
                j = find(i == interfacedata.evalVariables);
                variables = interfacedata.evalMap{j}.variableIndex;
            else
                variables = find(interfacedata.monomtable(i,:));
            end
            oJacobian(find(ismember(linear,variables))) = 1;
        end
%    end
end
if m > 0
    cJacobian = zeros(m,length(linear));
    for i = 1:size(F_struc,2)-1
        %if nonlinear(i)%variabletype(i) | ismembc(i,interfacedata.evalVariables)
            f = F_struc(:,i+1);
            for j = find(f(:)')%1:size(F_struc,1)
                %   if f(j)%F_struc(j,i+1)
                %                if variabletype(i) | ismembc(i,interfacedata.evalVariables)
                if ismembc(i,interfacedata.evalVariables)
                    k = find(i == interfacedata.evalVariables);
                    variables = interfacedata.evalMap{k}.variableIndex;
                elseif nonlinear(i)
                    variables = find(interfacedata.monomtable(i,:));
                else
                    variables = i;;
                end
                cJacobian(j,find(ismembc(linear,variables))) = 1;
            end
       % else

       % end
    end
    [Nbegcol,Nlencol,Nrowndx] = lindosparse(cJacobian);
else
    cJacobian = [];
    Nbegcol = [];
    Nrowndx = [];
    Nlencol = [];
    Nbegcol = [Nbegcol sum(Nlencol)];
end
oJacobian = double(oJacobian | any(interfacedata.Q(linear,linear),2));

Nobjndx = find(oJacobian) - 1;
Nobjcnt = length(Nobjndx);
if  isempty(Nobjndx)
    Nobjndx = [];
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
filename='CantileverBeam_Quadrilateral_Bilinear_Structured';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'holes';
cost = {'compliance';'perimeter'};
weights = [1 0.1];
constraint = {'volumeConstraint'};
optimizer = 'SLERP';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';
constraints_case = 'INEQUALITY';
line_search_initiator = 'INCREASING LAST STEP';

nsteps = 10;
Vfrac_final = 0.2;
Perimeter_target=1;
optimality_final =1e-3;
constr_final =1e-3;

Vfrac_initial = 0.5;
optimality_initial = 1e-2;
constr_initial = 1e-3;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan