text:
function clout = reparse_continguous(cl)
% Re-define regions in region object based on contiguous blobs
%
% :Usage:
% ::
%
%    clout = reparse_continguous(cl)
%
% ..
%    NEEDS SOME ADDITIONAL WORK/CHECKING
% ..

ivec = region2imagevec(cl);

clout = region(ivec, 'contiguous_regions', 'noverbose');

% ----------------------------------------------
% Add other fields
% ----------------------------------------------

clindx = ivec.volInfo.cluster;  % wh are in new clusters
nvox = length(clindx);

dat = cat(2, cl.all_data);

if size(dat, 2) == nvox

    % matching field; get dat for new clout
    for j = 1:length(clout)
        clout(j).all_data = dat(:, clindx == j);

        clout(j).dat = nanmean(clout(j).all_data', 1)';
    end

elseif isempty(dat)
    % do nothing
else
    disp('Warning!  all_dat field is wrong size.');

end

end

author: canlab

repo: CanlabCore


text:
function sphere_llt_grid_display ( r, pc, lat_num, long_num, node_num, node_xyz, ...
    line_num, line_data, filename )

%*****************************************************************************80
%
%% SPHERE_LLT_GRID_DISPLAY displays points and lines of an LLT grid on a sphere.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the radius of the sphere.
%
%    Input, real PC(1,3), the center of the sphere.
%
%    Input, integer LAT_NUM, LONG_NUM, the number of latitude and longitude
%    lines to draw.  The latitudes do not include the North and South
%    poles, which will be included automatically, so LAT_NUM = 5, for instance,
%    will result in points along 7 lines of latitude.
%
%    Input, integer NODE_NUM, the number of grid points.
%
%    Input, real NODE_XYZ(NODE_NUM,3), the grid points.
%
%    Input, integer LINE_NUM, the number of grid lines.
%
%    Input, integer LINE_DATA(LINE_NUM,2), contains pairs of point indices for
%    line segments that make up the grid.
%
%    Input, string FILENAME, the name of a file in which to store a copy of the plot.
%

%
%  Get the scale.
%
  xyz_min(1) = min ( node_xyz(:,1) );
  xyz_max(1) = max ( node_xyz(:,1) );

  xyz_min(2) = min ( node_xyz(:,2) );
  xyz_max(2) = max ( node_xyz(:,2) );

  xyz_min(3) = min ( node_xyz(:,3) );
  xyz_max(3) = max ( node_xyz(:,3) );

  xyz_range(1:3) = xyz_max(1:3) - xyz_min(1:3);

  margin = 0.025 * max ( xyz_range(1), ...
                   max ( xyz_range(2), xyz_range(3) ) );

  x_min = xyz_min(1) - margin;
  x_max = xyz_max(1) + margin;
  y_min = xyz_min(2) - margin;
  y_max = xyz_max(2) + margin;
  z_min = xyz_min(3) - margin;
  z_max = xyz_max(3) + margin;
%
%  Draw the picture.
%
  figure ( )
  clf
  hold on
  point_size = 50;
  point_color = [ 0.0, 1.0, 0.0 ];
  scatter3 ( node_xyz(:,1), node_xyz(:,2), node_xyz(:,3), point_size, ...
    'k', 'filled' );

  for i = 1 : line_num
    p = [ line_data(i,1), line_data(i,2) ];
    line ( node_xyz(p,1), node_xyz(p,2), node_xyz(p,3), 'LineWidth', 2 );
  end
%
%  We want to include a sphere in the picture.  But if the sphere is close to
%  the correct radius, lines that should lie on the surface because they are
%  curved will actually tunnel through the surface because we draw them straight.
%
  [ x, y, z ] = sphere ( 40 );
  x = 0.95 * r * x;
  y = 0.95 * r * y;
  z = 0.95 * r * z;

  c = z ./ z;
  surf ( x, y, z, c, 'EdgeColor', 'None' );

  axis equal
  grid on
  xlabel ( '--X axis--' )
  ylabel ( '--Y axis--' )
  zlabel ( '--Z axis--' )
  title ( filename, 'FontSize', 24 )
  view ( 3 )
  hold off

  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saving plot as "%s".\n', filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function CI = calcCI(resp,time,censoring)
% resp: Hazards ratio for Cox regression
% Y: Time-to-event (the higher, the better)

% FLIP THE RESPONSE. If Resp1 > resp2 && Y1 > Y2, add + 1. IF Y1<Y2 and Y1
% is censored, discard.
% ________________________________________________________________________

% INITIALIZATION
nInst = numel(time);
resp = - resp; % Flipping the hazards for easier comparison

% GETTING ALL PAIRS OF RESPONSES
pairs = combnk(1:nInst,2);
resp = resp(pairs);
time = time(pairs);
censoring = censoring(pairs);

% REMOVING UNUSABLE DATA

% 1. Both patients lived: unusable data
sumCens = censoring(:,1) + censoring(:,2);
out = find(sumCens == 2);
resp(out,:) = []; time(out,:) = []; censoring(out,:) = [];

% 2. One patient died, but the follow-up of the other is not long enough (did not outlive the first patient): unusable data
valCheck = (censoring(:,2) - censoring(:,1)) .* (time(:,2) - time(:,1));
out = find(valCheck < 0);
resp(out,:) = []; time(out,:) = []; censoring(out,:) = [];

% 3. Both patients died at the same time: unusable data
valCheck1 = time(:,2) - time(:,1); valCheck2 = censoring(:,2) - censoring(:,1); % (all 1,1 censored pairs have been removed by now, leaving only 0,1 or 0,0)
out = find(valCheck1 == 0 & valCheck2 == 0);
resp(out,:) = []; time(out,:) = []; censoring(out,:) = [];

% COMPUTING TIES IN RESPONSES
diffResp = resp(:,2) - resp(:,1);
indTies = find(diffResp == 0);
nTies = numel(indTies);
resp(indTies,:) = []; time(indTies,:) = []; censoring(indTies,:) = [];

% CALCULATING CI
valCheck = (time(:,2) - time(:,1)) .* (resp(:,2) - resp(:,1)); nCheck = numel(valCheck);
indGood = find(valCheck > 0); nGood = numel(indGood);
CI = (nGood + 0.5*nTies)/(nCheck + nTies);
if isnan(CI)
    CI = 0.5;
end

end

author: mvallieres

repo: radiomics


text:
%% Test constant objective

%% LP
clc
% Linear Objective & Bias
f = -[6 5]';
objbias = 5;
% Linear Constraints
A = ([1,4; 6,4; 2, -5]);
b = [16;28;6];
lb = [0;0]; ub = [10;10];
%Build Object
Opt = opti('f',f,'objbias',objbias,'ineq',A,b,'bounds',lb,ub)
%Build & Solve
[x,fval,exitflag,info] = solve(Opt)
%Plot
plot(Opt)

%% QP
P = coinRead('testQP.qps')

O = opti(P,optiset('solver','scip','display','iter'))
[x,f] = solve(O)
plot(O)

%% MILP2
clc
%Objective & Constraints
f = -[1 2 3 1]';
A = [-1 1 1 10; 1 -3 1 0];
b = [20;30];
Aeq = [0 1 0 -3.5];
beq = 0;
%Setup Options
opts = optiset('solver','cbc','display','iter');
%Build & Solve
Opt = opti('grad',f,'ineq',A,b,'eq',Aeq,beq,'bounds',[0 0 0 2]',[40 inf inf 3]','int','CCCI','options',opts,'objbias',10)
[x,fval,exitflag,info] = solve(Opt)
plot(Opt)

%% QCQP 1
clc
%Objective & Constraints
H = [33 6    0;
     6  22   11.5;
     0  11.5 11];
f = [-1;-2;-3];
A = [-1 1 1; 1 -3 1];
b = [20;30];
Q = eye(3);
l = [2;2;2];
r = 0.5;
lb = [0;0;0];
ub = [40;inf;inf];
%Build & Solve
Opt = opti('H',H,'f',f,'ineq',A,b,'qc',Q,l,r,'bounds',lb,ub,'objbias',10,'options',optiset('display','iter','solver','ipopt'))
[x,fval,exitflag,info] = solve(Opt)

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function tests = test_readMesh
% Test suite for the file readMesh.
%
%   Test suite for the file readMesh
%
%   Example
%   test_readMesh
%
%   See also
%     readMesh

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2021-01-25,    using Matlab 9.8.0.1323502 (R2020a)
% Copyright 2021 INRAE - BIA-BIBS.

tests = functiontests(localfunctions);

function test_off_VF(testCase) %#ok<*DEFNU>

[v, f] = readMesh('mushroom.off');

assertEqual(testCase, 226, size(v, 1));
assertEqual(testCase, 448, size(f, 1));

function test_off_single(testCase) %#ok<*DEFNU>

mesh = readMesh('mushroom.off');

assertTrue(testCase, isstruct(mesh));
assertTrue(testCase, isfield(mesh, 'vertices'));
assertTrue(testCase, isfield(mesh, 'faces'));

assertEqual(testCase, 226, size(mesh.vertices, 1));
assertEqual(testCase, 448, size(mesh.faces, 1));

function test_ply_VF(testCase) %#ok<*DEFNU>

[v, f] = readMesh('apple.ply');

assertEqual(testCase, 867, size(v, 1));
assertEqual(testCase, 1704, size(f, 1));

function test_ply_single(testCase) %#ok<*DEFNU>

mesh = readMesh('apple.ply');

assertTrue(testCase, isstruct(mesh));
assertTrue(testCase, isfield(mesh, 'vertices'));
assertTrue(testCase, isfield(mesh, 'faces'));

assertEqual(testCase, 867, size(mesh.vertices, 1));
assertEqual(testCase, 1704, size(mesh.faces, 1));

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [A,b,x] = foxgood(n)
%FOXGOOD Test problem: severely ill-posed problem.
%
% [A,b,x] = foxgood(n)
%
% This is a model problem which does not satisfy the
% discrete Picard condition for the small singular values.
% The problem was first used by Fox & Goodwin.

% Reference: C. T. H. Baker, "The Numerical Treatment of
% Integral Equations", Clarendon Press, Oxford, 1977; p. 665.

% Discretized by simple quadrature (midpoint rule).

% Per Christian Hansen, IMM, 03/16/93.

% Initialization.
h = 1/n; t = h*((1:n)' - 0.5);

A = h*sqrt((t.^2)*ones(1,n) + ones(n,1)*(t.^2)');
x = t; b = ((1+t.^2).^1.5 - t.^3)/3;

repo: REST


text:
function D = spm_eeg_inv_vbecd_gui(D,val)
% GUI function for Bayesian ECD inversion
% - load the necessary data, if not provided
% - fill in all the necessary bits for the VB-ECD inversion routine,
% - launch the B_ECD routine, aka. spm_eeg_inv_vbecd
% - displays the results.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
%
% $Id: spm_eeg_inv_vbecd_gui.m 7702 2019-11-22 11:32:26Z guillaume $

%%
% Load data, if necessary
%==========
if nargin<1
    D = spm_eeg_load;
end

DISTTHRESH=120; %% mm distance from centre of brain classified as inside the head

%%
% Check if the forward model was prepared & handle the other info bits
%========================================
if ~isfield(D,'inv')
    error('Data must have been prepared for inversion procedure...')
end
if nargin==2
    % check index provided
    if val>length(D.inv)
        val = length(D.inv);
        D.val = val;
    end
else
    if isfield(D,'val')
        val = D.val;
    else
        % use last one
        val = length(D.inv);
        D.val = val;
    end
end

% Use val to define which is the "current" inv{} to use
% If no inverse solution already calculated (field 'inverse' doesn't exist)
% use that inv{}. Otherwise create a new one by copying the previous
% inv{} structure
if isfield(D.inv{val},'inverse')
    % create an extra inv{}
    Ninv = length(D.inv);
    D.inv{Ninv+1} = D.inv{val};
    if isfield(D.inv{Ninv+1},'contrast')
        % no contrast field used here !
        D.inv{Ninv+1} = rmfield(D.inv{Ninv+1},'contrast');
    end
    val = Ninv+1;
    D.val = val;
end

if ~isfield(D.inv{val}, 'date')
    % Set time , date, comments & modality
    clck = fix(clock);
    if clck(5) < 10
        clck = [num2str(clck(4)) ':0' num2str(clck(5))];
    else
        clck = [num2str(clck(4)) ':' num2str(clck(5))];
    end
    D.inv{val}.date = strvcat(date,clck);
end

if ~isfield(D.inv{val}, 'comment'),
   D.inv{val}.comment = {spm_input('Comment/Label for this analysis', '+1', 's')};
end

D.inv{val}.method = 'vbecd';

%% Struct that collects the inputs for vbecd code
P = [];

P.modality = spm_eeg_modality_ui(D, 1, 1);

for m = 1:numel(D.inv{val}.datareg)
    if strncmp(D.inv{val}.datareg(m).modality, P.modality, 3)
        modind = m;
    end
end

data = spm_eeg_inv_get_vol_sens(D, val, [], 'inv', P.modality);

P.forward.vol     = data.(P.modality(1:3)).vol;
if ischar(P.forward.vol)
    P.forward.vol = ft_read_headmodel(P.forward.vol);
end

P.forward.sens    = data.(P.modality(1:3)).sens;
P.forward.siunits = data.siunits;
M1 = data.transforms.toMNI;
[U, L, V] = svd(M1(1:3, 1:3));
orM1(1:3,1:3) =U*V'; %% for switching orientation between meg and mni space

P.Ic = indchantype(D, P.modality, 'GOOD');

if isempty(P.Ic)
    error(['The specified modality (' P.modality ') is missing from file ' D.fname]);
else
    P.channels = D.chanlabels(P.Ic);
end

P.forward.chanunits = D.units(P.Ic);

[P.forward.vol, P.forward.sens] =  ft_prepare_vol_sens( ...
    P.forward.vol, P.forward.sens, 'channel', P.channels);

P.forward.sens.prj = D.coor2D(P.Ic);

%%
% Deal with data
%===============

% time bin or time window
msg_tb = ['time_bin or average_win [',num2str(round(min(D.time)*1e3)), ...
            ' ',num2str(round(max(D.time)*1e3)),'] ms'];
ask_tb = 1;
while ask_tb
    tb = spm_input(msg_tb,1,'r');   % ! in msec
    if length(tb)==1
        if tb>=min(D.time([], 'ms')) && tb<=max(D.time([], 'ms'))
            ask_tb = 0;
        end
    elseif length(tb)==2
        if all(tb>=floor(min(D.time([], 'ms')))) && all(tb<=ceil(max(D.time([], 'ms')))) && tb(1)<=tb(2)
            ask_tb = 0;
        end
    end
end
if length(tb)==1
    [kk,ltb] = min(abs(D.time([], 'ms')-tb)); % round to nearest time bin
else
    [kk,ltb(1)] = min(abs(D.time([], 'ms')-tb(1)));  % round to nearest time bin
    [kk,ltb(2)] = min(abs(D.time([], 'ms')-tb(2)));
    ltb = ltb(1):ltb(2); % list of time bins 'tb' to use
end

% trial type
if D.ntrials>1
    msg_tr = ['Trial type number [1 ',num2str(D.ntrials),']'];
    ltr = spm_input(msg_tr,2,'i',num2str(1:D.ntrials));
    tr_q = 1;
else
    tr_q = 0;
    ltr = 1;
end

% data, averaged over time window considered

EEGscale=1;

%% SORT OUT EEG UNITS AND CONVERT VALUES TO VOLTS
if strcmpi(P.modality,'EEG')
    allunits  = strvcat('uV','mV','V');
    allscales = [1, 1e3, 1e6]; %%
    EEGscale  = 0;
    eegunits  = unique(D.units(D.indchantype('EEG')));
    Neegchans = numel(D.units(D.indchantype('EEG')));
    for j=1:length(allunits)
        if strcmp(deblank(allunits(j,:)),deblank(eegunits))
            EEGscale=allscales(j);
        end % if
    end % for j

    if EEGscale==0
        warning('units unspecified');
        if mean(std(D(P.Ic,ltb,ltr)))>1e-2
            guess_ind=[1 2 3];
        else
            guess_ind=[3 2 1];
        end
        msg_str = sprintf('Units of EEG are %s ? (rms=%3.2e)',allunits(guess_ind(1),:),mean(std(D(P.Ic,ltb,ltr))));
        dip_ch  = sprintf('%s|%s|%s',allunits(guess_ind(1),:),allunits(guess_ind(2),:),allunits(guess_ind(3),:));
        dip_val = [1,2,3];
        def_opt = 1;
        unitind = spm_input(msg_str,2,'b',dip_ch,dip_val,def_opt);
        %ans=spm_input(msg_str,1,'s','yes');
        allunits(guess_ind(unitind),:)
        D = units(D, 1:Neegchans, allunits(guess_ind(unitind),:));
        EEGscale = allscales(guess_ind(unitind));
        D.save; %% Save the new units
    end %if EEGscale==0

end % if eeg data

dat_y = squeeze(mean(D(P.Ic,ltb,ltr)*EEGscale,2));

%%
% Other bits of the P structure, apart for priors and #dipoles
%==============================

P.ltr          = ltr;
P.Nc           = length(P.Ic);

% Deal with dipoles number and priors
%====================================
dip_q = 0; % number of dipole 'elements' added (single or pair)
dip_c = 0; % total number of dipoles in the model
adding_dips = 1;
clear dip_pr

priorlocvardefault=[100, 100, 100]; %% location variance default in mm
nopriorlocvardefault=[80*80, 80*80, 80*80];
nopriormomvardefault=[10, 10, 10]*100; %% moment variance in nAM
priormomvardefault=[1, 1, 1]; %%

while adding_dips
    if dip_q>0,
        msg_dip =['Add dipoles to ',num2str(dip_c),' or stop?'];
        dip_ch = 'Single|Symmetric Pair|Stop';
        dip_val = [1,2,0];
        def_opt=3;
    else
        msg_dip =['Add dipoles to model'];
        def_opt=1;
        dip_ch = 'Single|Symmetric Pair';
        dip_val = [1,2];
    end
    a_dip = spm_input(msg_dip,2+tr_q+dip_q,'b',dip_ch,dip_val,def_opt);
    if a_dip == 0
        adding_dips = 0;
    elseif a_dip == 1
    % add a single dipole to the model
        dip_q = dip_q+1;
        dip_pr(dip_q) = struct( 'a_dip',a_dip, ...
            'mu_w0',[],'mu_s0',[],'S_s0',eye(3),'S_w0',eye(3));
          %%  'ab20',[],'ab30',[]); %% ,'Tw', [],'Ts', []);
        % Location prior
        spr_q = spm_input('Location prior ?',1+tr_q+dip_q+1,'b', ...
                    'Informative|Non-info',[1,0],2);
        if spr_q
            % informative location prior
            str = 'Location prior';
            while 1
                s0mni = spm_input(str, 1+tr_q+dip_q+2,'e',[0 0 0])';

                %outside = ~ft_inside_headmodel(s0mni',mnivol);

                distcentre=sqrt(dot(s0mni'-mean(D.inv{D.val}.mesh.tess_mni.vert),s0mni'-mean(D.inv{D.val}.mesh.tess_mni.vert))); %%
                outside=(distcentre>DISTTHRESH);

                s0=D.inv{val}.datareg(modind).fromMNI*[s0mni' 1]';
                s0=s0(1:3);

                str2='Prior location variance (mm2)';
                diags_s0_mni = spm_input(str2, 1+tr_q+dip_q+2,'e',priorlocvardefault)';

                S_s0_ctf=orM1*diag(diags_s0_mni)*orM1'; %% transform covariance

                %% need to leave diags(S0) free

                if all(~outside), break, end
                    str = 'Prior location must be inside head';
                  end
            dip_pr(dip_q).mu_s0 = s0;
           else
            % no location  prior
            dip_pr(dip_q).mu_s0 = zeros(3,1);
            diags_s0_mni= nopriorlocvardefault';
            S_s0_ctf=diag(diags_s0_mni);
            end

        dip_pr(dip_q).S_s0=S_s0_ctf; %

        % Moment prior
        wpr_q = spm_input('Moment prior ?',1+tr_q+dip_q+spr_q+2,'b', ...
                    'Informative|Non-info',[1,0],2);
        if wpr_q
            % informative moment prior
            w0_mni= spm_input('Moment prior', ...
                                        1+tr_q+dip_q+spr_q+3,'e',[0 0 0])';
            str2='Prior moment variance (nAm2)';
            diags_w0_mni = spm_input(str2, 1+tr_q+dip_q+2,'e',priormomvardefault)';
            dip_pr(dip_q).mu_w0 =orM1*w0_mni;
            S_w0_ctf=orM1*diag(diags_w0_mni)*orM1';

        else
            % no location  prior
            dip_pr(dip_q).mu_w0 = zeros(3,1);
            S_w0_ctf= diag(nopriormomvardefault);
        end
        %% set up covariance matrix for orientation with no crosstalk terms (for single
        %% dip)
        dip_pr(dip_q).S_w0=S_w0_ctf;
        dip_c = dip_c+1;
    else
    % add a pair of symmetric dipoles to the model
        dip_q = dip_q+1;
        dip_pr(dip_q) = struct( 'a_dip',a_dip, ...
            'mu_w0',[],'mu_s0',[],'S_s0',eye(6),'S_w0',eye(6));
        %%...
          %  'ab20',[],'ab30',[]); %%,'Tw',eye(6),'Ts',eye(6));
        % Location prior
        spr_q = spm_input('Location prior ?',1+tr_q+dip_q+1,'b', ...
                    'Informative|Non-info',[1,0],2);
        if spr_q
            % informative location prior
             str = 'Location prior (one side only)';
            while 1
                s0mni = spm_input(str, 1+tr_q+dip_q+2,'e',[0 0 0])';
                syms0mni=s0mni;
                syms0mni(1)=-syms0mni(1);

                distcentre=sqrt(dot(s0mni'-mean(D.inv{D.val}.mesh.tess_mni.vert),s0mni'-mean(D.inv{D.val}.mesh.tess_mni.vert))); %%
                outside=(distcentre>DISTTHRESH);

                s0=D.inv{val}.datareg(modind).fromMNI*[s0mni' 1]';

                s0sym=D.inv{val}.datareg(modind).fromMNI*[syms0mni' 1]';

                str2='Prior location variance (mm2)';
                tmp_diags_s0_mni = spm_input(str2, 1+tr_q+dip_q+2,'e',priorlocvardefault)';
                tmp_diags_s0_mni= [tmp_diags_s0_mni ; tmp_diags_s0_mni ];

                if all(~outside), break, end
                    str = 'Prior location must be inside head';
                  end

            dip_pr(dip_q).mu_s0 = [s0(1:3);s0sym(1:3)];
        else
            % no location  prior
            dip_pr(dip_q).mu_s0 = zeros(6,1);
            tmp_diags_s0 = [nopriorlocvardefault';nopriorlocvardefault'];
            tmp_diags_s0_mni=[nopriorlocvardefault';nopriorlocvardefault'];
        end %% end of if informative prior
        %% setting up a covariance matrix where there is covariance between
        %% the x parameters negatively coupled, y,z positively.
         mni_dip_pr(dip_q).S_s0 = eye(length(tmp_diags_s0_mni)).*repmat(tmp_diags_s0_mni,1,length(tmp_diags_s0_mni));
         mni_dip_pr(dip_q).S_s0(4,1)=-mni_dip_pr(dip_q).S_s0(4,4); % reflect in x
         mni_dip_pr(dip_q).S_s0(5,2)=mni_dip_pr(dip_q).S_s0(5,5); % maintain y and z
         mni_dip_pr(dip_q).S_s0(6,3)=mni_dip_pr(dip_q).S_s0(6,6);

         mni_dip_pr(dip_q).S_s0(1,4)=mni_dip_pr(dip_q).S_s0(4,1);
         mni_dip_pr(dip_q).S_s0(2,5)=mni_dip_pr(dip_q).S_s0(5,2);
         mni_dip_pr(dip_q).S_s0(3,6)=mni_dip_pr(dip_q).S_s0(6,3);
         %% transform to MEG space

         %dip_pr(dip_q).S_s0(:,1:3)=orM1*mni_dip_pr(dip_q).S_s0(:,1:3)*orM1'; %% NEED TO LOOK AT THIS
         %dip_pr(dip_q).S_s0(:,4:6)=orM1*mni_dip_pr(dip_q).S_s0(:,4:6)*orM1';
         tmp1=orM1*mni_dip_pr(dip_q).S_s0(1:3,1:3)*orM1'; %% NEED TO LOOK AT THIS
         tmp2=orM1*mni_dip_pr(dip_q).S_s0(1:3,4:6)*orM1';
         tmp3=orM1*mni_dip_pr(dip_q).S_s0(4:6,4:6)*orM1';
         tmp4=orM1*mni_dip_pr(dip_q).S_s0(4:6,1:3)*orM1';
         dip_pr(dip_q).S_s0(1:3,1:3)=tmp1;
         dip_pr(dip_q).S_s0(1:3,4:6)=tmp2;
         dip_pr(dip_q).S_s0(4:6,4:6)=tmp3;
         dip_pr(dip_q).S_s0(4:6,1:3)=tmp4;

        % Moment prior
        wpr_q = spm_input('Moment prior ?',1+tr_q+dip_q+spr_q+2,'b', ...
                                           'Informative|Non-info',[1,0],2);
        if wpr_q
            % informative moment prior
            tmp= spm_input('Moment prior (right only)', ...
                                      1+tr_q+dip_q+spr_q+3,'e',[1 1 1])';
            tmp = [tmp ; tmp] ; tmp(4) = tmp(4);
            dip_pr(dip_q).mu_w0 = tmp;

            str2='Prior moment variance (nAm2)';
            diags_w0 = spm_input(str2, 1+tr_q+dip_q+spr_q+3,'e',priormomvardefault)';
            tmp_diags_w0=[diags_w0; diags_w0];

        else
            % no moment  prior
            dip_pr(dip_q).mu_w0 = zeros(6,1);
            tmp_diags_w0 = [nopriormomvardefault'; nopriormomvardefault'];
        end
        %dip_pr(dip_q).S_w0=eye(length(diags_w0)).*repmat(diags_w0,1,length(diags_w0));
        %% couple all orientations, except x, positively or leave for now...

                mni_dip_pr(dip_q).S_w0 = eye(length(tmp_diags_w0)).*repmat(tmp_diags_w0,1,length(tmp_diags_w0));

            mni_dip_pr(dip_q).S_w0(4,1)=-mni_dip_pr(dip_q).S_w0(4,4); % reflect x orientation
            mni_dip_pr(dip_q).S_w0(5,2)=mni_dip_pr(dip_q).S_w0(5,5); %
            mni_dip_pr(dip_q).S_w0(6,3)=mni_dip_pr(dip_q).S_w0(6,6); %
            mni_dip_pr(dip_q).S_w0(1,4)=-mni_dip_pr(dip_q).S_w0(4,1); %
            mni_dip_pr(dip_q).S_w0(2,5)=mni_dip_pr(dip_q).S_w0(5,2); %
            mni_dip_pr(dip_q).S_w0(3,6)=mni_dip_pr(dip_q).S_w0(6,3); %
            tmp1=orM1*mni_dip_pr(dip_q).S_w0(1:3,1:3)*orM1';
            tmp2=orM1*mni_dip_pr(dip_q).S_w0(1:3,4:6)*orM1';
            tmp3=orM1*mni_dip_pr(dip_q).S_w0(4:6,4:6)*orM1';
            tmp4=orM1*mni_dip_pr(dip_q).S_w0(4:6,1:3)*orM1';
            dip_pr(dip_q).S_w0(1:3,1:3)=tmp1;
            dip_pr(dip_q).S_w0(1:3,4:6)=tmp2;
            dip_pr(dip_q).S_w0(4:6,4:6)=tmp3;
            dip_pr(dip_q).S_w0(4:6,1:3)=tmp4;

        dip_c = dip_c+2;
    end
end

%str2='Data SNR (amp)';
% SNRamp = spm_input(str2, 1+tr_q+dip_q+2+1,'e',5)';
%SNRamp=3;
SNRamp=3;
 hE=log(SNRamp^2); %% expected log precision of data
 hC=1; % variability of the above precision

str2='Number of iterations';
 Niter = spm_input(str2, 1+tr_q+dip_q+2+2,'e',10)';

%%
% Get all the priors together and build structure to pass to inv_becd

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

priors = struct('mu_w0',cat(1,dip_pr(:).mu_w0), ...
                'mu_s0',cat(1,dip_pr(:).mu_s0), ...
                'S_w0',blkdiag(dip_pr(:).S_w0),'S_s0',blkdiag(dip_pr(:).S_s0),'hE',hE,'hC',hC);

P.priors = priors;

%%
% Launch inversion !
%===================

% Initialise inverse field
inverse = struct( ...
    'F',[], ... % free energy
    'pst',D.time, ... % all time points in data epoch
    'tb',tb, ... % time window/bin used
    'ltb',ltb, ... % list of time points used
    'ltr',ltr, ... % list of trial types used
    'n_seeds',length(ltr), ... % using this field for multiple reconstruction
    'n_dip',dip_c, ... % number of dipoles used
    'loc',[], ... % loc of dip (3 x n_dip)
    'j',[], ... % dipole(s) orient/ampl, in 1 column
    'cov_loc',[], ... % cov matrix of source location
    'cov_j',[], ...   % cov matrix of source orient/ampl
    'Mtb',1, ... % ind of max EEG power in time series, 1 as only 1 tb.
    'exitflag',[], ... % Converged (1) or not (0)
    'P',[]);             % save all kaboodle too.

for ii=1:length(ltr)
    P.y = dat_y(:,ii);
    P.ii = ii;

 %% set up figures
    P.handles.hfig  = spm_figure('GetWin','Graphics');
    spm_clf(P.handles.hfig)
    P.handles.SPMdefaults.col = get(P.handles.hfig,'colormap');
    P.handles.SPMdefaults.renderer = get(P.handles.hfig,'renderer');
    set(P.handles.hfig,'userdata',P)
    dip_amp=[];
    for j=1:Niter,

     Pout(j) = spm_eeg_inv_vbecd(P);
     close(gcf);
     varresids(j)=var(Pout(j).y-Pout(j).ypost);
     pov(j)=100*(1-varresids(j)/var(Pout(j).y)); %% percent variance explained
     allF(j)=Pout(j).F;
     dip_mom=reshape(Pout(j).post_mu_w,3,length(Pout(j).post_mu_w)/3);
     dip_amp(j,:)=sqrt(dot(dip_mom,dip_mom));
    % display
     megloc=reshape(Pout(j).post_mu_s,3,length(Pout(j).post_mu_s)/3); % loc of dip (3 x n_dip)
     mniloc=D.inv{val}.datareg(modind).toMNI*[megloc;ones(1,size(megloc,2))]; %% actual MNI location (with scaling)
     megmom=reshape(Pout(j).post_mu_w,3,length(Pout(j).post_mu_w)/3); % moments of dip (3 x n_dip)
     megposvar=reshape(diag(Pout(j).post_S_s),3,length(Pout(j).post_mu_s)/3); %% estimate of positional uncertainty in three principal axes
     mnimom=orM1*megmom; %% convert moments into mni coordinates through a rotation (no scaling or translation)
     mniposvar=(orM1*sqrt(megposvar)).^2; %% convert pos variance into approx mni space by switching axes

     displayVBupdate2(Pout(j).y,pov,allF,Niter,dip_amp,mnimom,mniloc(1:3,:),mniposvar,P,j,[],Pout(j).F,Pout(j).ypost,[]);

    end; % for  j
    allF=[Pout.F];
    [maxFvals,maxind]=max(allF);

    P=Pout(maxind); %% take best F
    % Get the results out.
    inverse.pst = tb*1e3;
    inverse.F(ii) = P.F; % free energy

     megloc=reshape(P.post_mu_s,3,length(P.post_mu_s)/3); % loc of dip (3 x n_dip)
     meg_w=reshape(P.post_mu_w,3,length(P.post_mu_w)/3); % moments of dip (3 x n_dip)
     mni_w=orM1*meg_w; %% orientation in mni space
     mniloc=D.inv{val}.datareg(modind).toMNI*[megloc;ones(1,size(megloc,2))]; %% actual MNI location (with scaling)
     inverse.mniloc{ii}=mniloc(1:3,:);
    inverse.loc{ii} = megloc;

    inverse.j{ii} = P.post_mu_w; % dipole(s) orient/ampl, in 1 column in meg space
    inverse.jmni{ii} = reshape(mni_w,1,prod(size(mni_w)))'; % dipole(s) orient/ampl in mni space
    inverse.cov_loc{ii} = P.post_S_s; % cov matrix of source location
    inverse.cov_j{ii} = P.post_S_w; % cov matrix of source orient/ampl
    inverse.exitflag(ii) = 1; % Converged (1) or not (0)
    inverse.P{ii} = P; % save all kaboodle too.
    %% show final result
    pause(1);

    spm_clf(P.handles.hfig)

     megmom=reshape(Pout(maxind).post_mu_w,3,length(Pout(maxind).post_mu_w)/3); % moments of dip (3 x n_dip)
     %megposvar=reshape(diag(Pout(maxind).post_S_s),3,length(Pout(maxind).post_mu_s)/3); %% estimate of positional uncertainty in three principal axes
     mnimom=orM1*megmom; %% convert moments into mni coordinates through a rotation (no scaling or translation)

     longorM1=zeros(size(Pout(maxind).post_S_s,1));
     for k1=1:length(Pout(maxind).post_S_s)/3;
        longorM1((k1-1)*3+1:k1*3,(k1-1)*3+1:k1*3)=orM1;
     end; % for k1
     S0_mni=longorM1*Pout(maxind).post_S_s*longorM1';
     mniposvar=diag(S0_mni); %% convert pos variance into approx mni space by switching axes
     mniposvar=reshape(mniposvar,3,length(Pout(maxind).post_S_s)/3);

     displayVBupdate2(Pout(maxind).y,pov,allF,Niter,dip_amp,mnimom,mniloc(1:3,:),mniposvar,P,j,[],Pout(maxind).F,Pout(maxind).ypost,maxind);

    %displayVBupdate2(Pout(maxind).y,pov,allF,Niter,dip_amp,mniloc,Pout(maxind).post_mu_s,Pout(maxind).post_S_s,P,j,[],Pout(maxind).F,Pout(maxind).ypost,maxind,D);
  %
end

%  [tmp] = ft_compute_leadfield([-42 15 48], P.forward.sens, P.forward.vol,'reducerank',2);
%  estfield=tmp*[-21 -8 -45]';
%  plot(P.y,estfield./1000,'x')
%  axis equal

D.inv{val}.inverse = inverse;

%%
% Save results and display
%-------------------------
save(D)

return

function [P] = displayVBupdate2(y,pov_iter,F_iter,maxit,dipamp_iter,mu_w,mu_s,diagS_s,P,it,flag,F,yHat,maxind)

%% yHat is estimate of y based on dipole position
if ~exist('flag','var')
    flag = [];
end
if ~exist('maxind','var')
    maxind = [];
end

if isempty(flag) || isequal(flag,'ecd')
    % plot dipoles
    try
        opt.ParentAxes = P.handles.axesECD;
        opt.hfig = P.handles.hfig;
        opt.handles.hp = P.handles.hp;
        opt.handles.hq = P.handles.hq;
        opt.handles.hs = P.handles.hs;
        opt.handles.ht = P.handles.ht;
        opt.query = 'replace';
    catch
        P.handles.axesECD = axes(...
            'parent',P.handles.hfig,...
            'Position',[0.13 0.55 0.775 0.4],...
            'hittest','off',...
            'visible','off',...
            'deleteFcn',@back2defaults);
        opt.ParentAxes = P.handles.axesECD;
        opt.hfig = P.handles.hfig;
    end
    w = reshape(mu_w,3,[]);
    s = reshape(mu_s, 3, []);
%     mesh.faces=D.inv{D.val}.forward.mesh.face; %% in ctf space
%     mesh.vertices=D.inv{D.val}.forward.mesh.vert; %% in ctf space
%     [out] = spm_eeg_displayECD_ctf(...
%         s,w,reshape(diag(S_s),3,[]),[],mesh,opt);
       [out] = spm_eeg_displayECD(...
         s,w,diagS_s,[],opt);

        P.handles.hp = out.handles.hp;
        P.handles.hq = out.handles.hq;
        P.handles.hs = out.handles.hs;
        P.handles.ht = out.handles.ht;

end

% plot data and predicted data
pos = P.forward.sens.prj;
ChanLabel = P.channels;
in.f = P.handles.hfig;
in.noButtons = 1;
try
    P.handles.axesY;
catch
    figure(P.handles.hfig)
    P.handles.axesY = axes(...
        'Position',[0.02 0.3 0.3 0.2],...
        'hittest','off');
    in.ParentAxes = P.handles.axesY;
    spm_eeg_plotScalpData(y,pos,ChanLabel,in);
    title(P.handles.axesY,'measured data')
end
if isempty(flag) || isequal(flag,'data') || isequal(flag,'ecd')
    %yHat = P.gmn*mu_w;
    miY = min([yHat;y]);
    maY = max([yHat;y]);
    try
        P.handles.axesYhat;
        d = get(P.handles.axesYhat,'userdata');
        yHat = yHat(d.goodChannels);
        clim = [min(yHat(:))-( max(yHat(:))-min(yHat(:)) )/63,...
            max(yHat(:))];
        ZI = griddata(...
            d.interp.pos(1,:),d.interp.pos(2,:),full(double(yHat)),...
            d.interp.XI,d.interp.YI);
        set(d.hi,'Cdata',flipud(ZI));
        caxis(P.handles.axesYhat,clim);
        delete(d.hc)
        [C,d.hc] = contour(P.handles.axesYhat,flipud(ZI),...
            'linecolor',0.5.*ones(3,1));
        set(P.handles.axesYhat,...
            'userdata',d);
    catch
        figure(P.handles.hfig)
        P.handles.axesYhat = axes(...
            'Position',[0.37 0.3 0.3 0.2],...
            'hittest','off');
        in.ParentAxes = P.handles.axesYhat;
        spm_eeg_plotScalpData(yHat,pos,ChanLabel,in);
        title(P.handles.axesYhat,'predicted data')
    end
    try
        P.handles.axesYhatY;
    catch
        figure(P.handles.hfig)
        P.handles.axesYhatY = axes(...
            'Position',[0.72 0.3 0.25 0.2],...
            'NextPlot','replace',...
            'box','on');
    end
    plot(P.handles.axesYhatY,y,yHat,'.')
    set(P.handles.axesYhatY,...
        'nextplot','add')
    plot(P.handles.axesYhatY,[miY;maY],[miY;maY],'r')
    set(P.handles.axesYhatY,...
        'nextplot','replace')
    title(P.handles.axesYhatY,'predicted vs measured data')
    axis(P.handles.axesYhatY,'square','tight')
    grid(P.handles.axesYhatY,'on')

end

if isempty(flag) || isequal(flag,'var')
    % plot precision hyperparameters
    try
        P.handles.axesVar1;
    catch
        figure(P.handles.hfig)
        P.handles.axesVar1 = axes(...
            'Position',[0.05 0.05 0.25 0.2],...
            'NextPlot','replace',...
            'box','on');
    end
    plot(P.handles.axesVar1,F_iter,'o-');
    if ~isempty(maxind),
        hold on;
        h=plot(P.handles.axesVar1,maxind,F_iter(maxind),'rd');
        set(h,'linewidth',4);
        end;
    set(P.handles.axesVar1,'Xlimmode','manual');
    set(P.handles.axesVar1,'Xlim',[1 maxit]);
    set(P.handles.axesVar1,'Xtick',1:maxit);
    set(P.handles.axesVar1,'Xticklabel',num2str([1:maxit]'));
    set(P.handles.axesVar1,'Yticklabel','');
    title(P.handles.axesVar1,'Free energy ')
    axis(P.handles.axesVar1,'square');
    set(P.handles.axesVar1,'Ylimmode','auto'); %,'tight')

    grid(P.handles.axesVar1,'on')

    try
        P.handles.axesVar2;
    catch
        figure(P.handles.hfig)
        P.handles.axesVar2 = axes(...
            'Position',[0.37 0.05 0.25 0.2],...
            'NextPlot','replace',...
            'box','on');
    end

    plot(P.handles.axesVar2,pov_iter,'*-')
    if ~isempty(maxind),
        hold on;
        h=plot(P.handles.axesVar2,maxind,pov_iter(maxind),'rd');
        set(h,'linewidth',4);
        end;
    set(P.handles.axesVar2,'Xlimmode','manual');
    set(P.handles.axesVar2,'Xlim',[1 maxit]);
    set(P.handles.axesVar2,'Xtick',1:maxit);
    set(P.handles.axesVar2,'Xticklabel',num2str([1:maxit]'));
    set(P.handles.axesVar2,'Ylimmode','manual'); %,'tight')
     set(P.handles.axesVar2,'Ylim',[0 100]);
    set(P.handles.axesVar2,'Ytick',[0:20:100]);
    set(P.handles.axesVar2,'Yticklabel',num2str([0:20:100]'));

    %set(P.handles.axesVar2,'Yticklabel','');
    title(P.handles.axesVar2,'Percent variance explained');
    axis(P.handles.axesVar2,'square');

    grid(P.handles.axesVar2,'on')

    try
        P.handles.axesVar3;
    catch
        figure(P.handles.hfig)
        P.handles.axesVar3 = axes(...
            'Position',[0.72 0.05 0.25 0.2],...
            'NextPlot','replace',...
            'box','on');
    end
    plot(P.handles.axesVar3,1:it,dipamp_iter','o-');
      if ~isempty(maxind),
        hold on;
        h=plot(P.handles.axesVar3,maxind,dipamp_iter(maxind,:)','rd');
        set(h,'linewidth',4);
        end;

    set(P.handles.axesVar3,'Xlimmode','manual');
    set(P.handles.axesVar3,'Xlim',[1 maxit]);
    set(P.handles.axesVar3,'Xtick',1:maxit);
    set(P.handles.axesVar3,'Xticklabel',num2str([1:maxit]'));
    set(P.handles.axesVar3,'Yticklabel','');
    title(P.handles.axesVar3,'Dipole amp (nAm) ')
    axis(P.handles.axesVar3,'square');
    set(P.handles.axesVar3,'Ylimmode','auto'); %,'tight')
    grid(P.handles.axesVar3,'on')

end

if ~isempty(flag) && (isequal(flag,'ecd') || isequal(flag,'mGN') )
    try
        P.handles.hte(2);
    catch
        figure(P.handles.hfig)
        P.handles.hte(2) = uicontrol('style','text',...
            'units','normalized',...
            'position',[0.2,0.91,0.6,0.02],...
            'backgroundcolor',[1,1,1]);
    end
    set(P.handles.hte(2),'string',...
        ['ECD locations: Modified Gauss-Newton scheme... ',num2str(floor(P.pc)),'%'])
else
    try
        set(P.handles.hte(2),'string','VB updates on hyperparameters')
    end
end

try
    P.handles.hte(1);
catch
    figure(P.handles.hfig)
    P.handles.hte(1) = uicontrol('style','text',...
        'units','normalized',...
        'position',[0.2,0.94,0.6,0.02],...
        'backgroundcolor',[1,1,1]);
end
try
    set(P.handles.hte(1),'string',...
        ['Model evidence: p(y|m) >= ',num2str(F(end),'%10.3e\n')])
end

try
    P.handles.hti;
catch
    figure(P.handles.hfig)
    P.handles.hti = uicontrol('style','text',...
        'units','normalized',...
        'position',[0.3,0.97,0.4,0.02],...
        'backgroundcolor',[1,1,1],...
        'string',['VB ECD inversion: trial #',num2str(P.ltr(P.ii))]);
end
drawnow

function back2defaults(e1,e2)
hf = spm_figure('FindWin','Graphics');
P = get(hf,'userdata');
try
    set(hf,'colormap',P.handles.SPMdefaults.col);
    set(hf,'renderer',P.handles.SPMdefaults.renderer);
end

author: spm


text:
function faces = compute_voronoi_triangulation(Q, vertex)

% compute_voronoi_triangulation - compute a triangulation
%
%   face = compute_voronoi_triangulation(Q);
%
%   Q is a Voronoi partition function, computed using
%       perform_fast_marching.
%   face(:,i) is the ith face.
%
%   Works in 2D and in 3D.
%
%   Copyright (c) 2008 Gabriel Peyre

d = nb_dims(Q);

if d==2
    faces = compute_voronoi_triangulation_2d(Q, vertex);
elseif d==3
    faces = compute_voronoi_triangulation_3d(Q, vertex);
else
    error('Works only for 2D or 3D data.');
end

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function faces = compute_voronoi_triangulation_3d(Q, vertex)

[dx dy dz] = meshgrid(0:1,0:1,0:1);
V = [];
for i=1:8
    v = Q(1+dx(i):end-1+dx(i),1+dy(i):end-1+dy(i),1+dz(i):end-1+dz(i));
    V = [V v(:)];
end

V = sort(V,2);
V = unique(V, 'rows');
V = V( prod(V,2)>0 ,:);

d = V(:,1)*0;
for i=1:7
   d = d+(V(:,i)~=V(:,i+1));
end
if sum(d==4)>1
%    warning('Problem with triangulation.');
end

% split 5 folds into 2
I = find(d==4);
faces = [];
for i=1:length(I)
    v = unique(V(I(i),:));
    x = vertex(:,v); % points
    % barycenter
    a = sum( x, 2 ) / 5;
    x = x-repmat(a, [1 size(x,2)]);
    t = atan2(x(2,:),x(1,:));
    [tmp,s] = sort(t);
    faces = [faces v( s([1 2 3 5]))'];
    faces = [faces v( s([3 4 1 5]))'];

end
%    faces = [V(I,1:3);V(I,[3 4 1])];

% remaining triangles
V = V( d==3 ,:);
for i=1:size(V,1)
    V(i,1:4) = unique(V(i,:));
end
faces = [faces, V(:,1:4)'];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function faces = compute_voronoi_triangulation_2d(Q, vertex)

V = [];
v = Q(1:end-1,1:end-1); V = [V v(:)];
v = Q(2:end,1:end-1); V = [V v(:)];
v = Q(1:end-1,2:end); V = [V v(:)];
v = Q(2:end,2:end); V = [V v(:)];

V = sort(V,2);
V = unique(V, 'rows');
V = V( prod(V,2)>0 ,:);

d = (V(:,1)~=V(:,2)) + (V(:,2)~=V(:,3)) + (V(:,3)~=V(:,4));
if sum(d==3)>1
    warning('Problem with triangulation.');
end

% split squares into 2 triangles
I = find(d==3);
faces = [];
for i=1:length(I)
    v = V(I(i),:);
    x = vertex(:,v); % points
    % barycenter
    a = sum( x, 2 ) / 4;
    x = x-repmat(a, [1 size(x,2)]);
    t = atan2(x(2,:),x(1,:));
    [tmp,s] = sort(t);
    faces = [faces; v( s([1 2 3]))];
    faces = [faces; v( s([3 4 1]))];

end
%    faces = [V(I,1:3);V(I,[3 4 1])];

% remaining triangles
V = V( d==2 ,:);
for i=1:size(V,1)
    V(i,1:3) = unique(V(i,:));
end
faces = [faces; V(:,1:3)];

% return in correct format
faces = faces';

author: alecjacobson

repo: gptoolbox


text:
function [handles] = playerturn(ha,handles)
%if it is the start of the turn then...
if handles.ipx==0
%     act only if the player is trying to move his own piece
    if handles.box(handles.x,handles.y)==handles.plr
        handles.ipx=handles.x;
        handles.ipy=handles.y;
    end
%if it is the completion of turn then...
else
    handles.fpx=handles.x;
    handles.fpy=handles.y;
    %first check whether its a valid turn or not.
    [check,handles.box] = checkfp(ha,handles.plr,handles.box,handles.ipx,handles.ipy,handles.fpx,handles.fpy,handles.dblbox,handles.blank);
    if check==1
%         if valid then update box
        handles.box(handles.ipx,handles.ipy)=0;
        handles.box(handles.fpx,handles.fpy)=handles.plr;
%         when someone's piece reaches its end line it is now double
        if handles.plr==1 && handles.fpx==8
            handles.dblbox(handles.fpx,handles.fpy)=1;
        end
        if handles.plr==2 && handles.fpx==1
            handles.dblbox(handles.fpx,handles.fpy)=1;
        end
%     if trying to move already doubled piece then update dblbox with new location
        if handles.dblbox(handles.ipx,handles.ipy)==1
            handles.dblbox(handles.ipx,handles.ipy)=0;
            handles.dblbox(handles.fpx,handles.fpy)=1;
        end
%         and change the display of the buttons accordingly
%         changing the initial button's string
        xy=handles.ipx*10+handles.ipy;
        set(ha(xy),'CData',handles.blank)
%         changing the destination button's string
        xy=handles.fpx*10+handles.fpy;
        if handles.dblbox(handles.fpx,handles.fpy)==1 && handles.plr==1
            set(ha(xy),'CData',handles.plr1dblmark)
        elseif handles.dblbox(handles.fpx,handles.fpy)==1 && handles.plr==2
            set(ha(xy),'CData',handles.plr2dblmark)
        elseif handles.plr==1
            set(ha(xy),'CData',handles.plr1mark)
        elseif handles.plr==2
            set(ha(xy),'CData',handles.plr2mark)
        end
%         and update player
        if handles.plr==1
            handles.plr=2;
            axes(ha(1));
            image(handles.plr2mark);
            axis off;
        else
            handles.plr=1;
            axes(ha(1));
            image(handles.plr1mark);
            axis off;

        end
%       now zero the initial position
        handles.ipx=0;
        handles.ipy=0;
%        loser command is applied after updating plr becoz one will turn
%        and other will lose, so ohter one needs to be checked
        loser=wholoses(handles.plr,handles.box);
        if loser~=0
            if loser==2
                msgbox('BLACK Wins','Winner','custom',handles.plr1mark);
            elseif loser==1
                msgbox('RED Wins','Winner','custom',handles.plr2mark);
            elseif loser==-1
                msgbox('Its a Draw');
            end
            close(handles.figure1);
        end
    else
%       if wrong input then initialize the positions so as to get input again
        msgbox('Wrong input','Error');
        handles.ipx=0;
        handles.ipy=0;
        handles.fpx=0;
        handles.fpy=0;
    end
end
end

author: Sable

repo: mcbench-benchmarks


text:
function com = pop_expevents(EEG, filename, unit)
% pop_expevents() - export events to CSV file
%
% Usage:
%   >> pop_expevents(EEG);
%   >> pop_expevents(EEG, filename);
%   >> pop_expevents(EEG, filename, unit);
%
% Inputs:
%   EEG      - EEGLAB dataset
%   filename - text file name
%   unit     - display latency and duration in 'samples' (default) or 'seconds'
%
% Outputs:
%   com      - Command to execute this function from the command line
%
% Examples:
%   Export all events and show a dialog window:
%     >> pop_expevents(EEG);
%
%   Export all events to a file, latency and duration in samples:
%     >> pop_expevents(EEG, filename);
%
%   Export all events to a file, latency and duration in seconds:
%     >> pop_expevents(EEG, filename, 'seconds');

% Copyright by Clemens Brunner <clbrunner@ucsd.edu>
% Revision: 0.10
% Date: 09/16/2011

% Revision history:
%   0.10: Initial version

% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 2 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 59 Temple
% Place - Suite 330, Boston, MA  02111-1307, USA.

com = '';
if nargin < 1
    help pop_expevents;
    return
end

if nargin < 2
	[filename, filepath] = uiputfile('*.csv', 'File name -- pop_expevents()');
    drawnow;
	if filename == 0
        return;  end
	filename = [filepath filename];
end;

if nargin < 3
    unit = 'samples';  end

if ~(strcmp(unit, 'samples') || strcmp(unit, 'seconds'))
    error('Unit must be either ''samples'' or ''seconds''.');  end

eeg_eventtable(EEG, 'dispTable', false, 'exportFile', filename, 'unit', unit);

com = sprintf('pop_expevents(%s, ''%s'', ''%s'');', inputname(1), filename, unit);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% normaliseiris - performs normalisation of the iris region by
% unwraping the circular region into a rectangular block of
% constant dimensions.
%
% Usage:
% [polar_array, polar_noise] = normaliseiris(image, x_iris, y_iris, r_iris,...
% x_pupil, y_pupil, r_pupil,eyeimage_filename, radpixels, angulardiv)
%
% Arguments:
% image                 - the input eye image to extract iris data from
% x_iris                - the x coordinate of the circle defining the iris
%                         boundary
% y_iris                - the y coordinate of the circle defining the iris
%                         boundary
% r_iris                - the radius of the circle defining the iris
%                         boundary
% x_pupil               - the x coordinate of the circle defining the pupil
%                         boundary
% y_pupil               - the y coordinate of the circle defining the pupil
%                         boundary
% r_pupil               - the radius of the circle defining the pupil
%                         boundary
% eyeimage_filename     - original filename of the input eye image
% radpixels             - radial resolution, defines vertical dimension of
%                         normalised representation
% angulardiv            - angular resolution, defines horizontal dimension
%                         of normalised representation
%
% Output:
% polar_array
% polar_noise
%
% Author:
% Libor Masek
% masekl01@csse.uwa.edu.au
% School of Computer Science & Software Engineering
% The University of Western Australia
% November 2003

function [polar_array, polar_noise] = normaliseiris(image, x_iris, y_iris, r_iris,...
x_pupil, y_pupil, r_pupil,eyeimage_filename, radpixels, angulardiv, method)

global DIAGPATH

radiuspixels = radpixels + 2;
angledivisions = angulardiv-1;

r = 0:(radiuspixels-1);

theta = 0:2*pi/angledivisions:2*pi;

x_iris = double(x_iris);
y_iris = double(y_iris);
r_iris = double(r_iris);

x_pupil = double(x_pupil);
y_pupil = double(y_pupil);
r_pupil = double(r_pupil);

% calculate displacement of pupil center from the iris center
ox = x_pupil - x_iris;
oy = y_pupil - y_iris;

if ox <= 0
    sgn = -1;
elseif ox > 0
    sgn = 1;
end

if ox==0 && oy > 0

    sgn = 1;

end

r = double(r);
theta = double(theta);

a = ones(1,angledivisions+1)* (ox^2 + oy^2);

% need to do something for ox = 0
if ox == 0
    phi = pi/2;
else
    phi = atan(oy/ox);
end

b = sgn.*cos(pi - phi - theta);

% calculate radius around the iris as a function of the angle
r = (sqrt(a).*b) + ( sqrt( a.*(b.^2) - (a - (r_iris^2))));

r = r - r_pupil;

rmat = ones(1,radiuspixels)'*r;

rmat = rmat.* (ones(angledivisions+1,1)*[0:1/(radiuspixels-1):1])';
rmat = rmat + r_pupil;

% exclude values at the boundary of the pupil iris border, and the iris scelra border
% as these may not correspond to areas in the iris region and will introduce noise.
%
% ie don't take the outside rings as iris data.
rmat  = rmat(2:(radiuspixels-1), :);

% calculate cartesian location of each data point around the circular iris
% region
xcosmat = ones(radiuspixels-2,1)*cos(theta);
xsinmat = ones(radiuspixels-2,1)*sin(theta);

xo = rmat.*xcosmat;
yo = rmat.*xsinmat;

xo = x_pupil+xo;
yo = y_pupil-yo;

% extract intensity values into the normalised polar representation through
% interpolation
[x,y] = meshgrid(1:size(image,2),1:size(image,1));
polar_array = interp2(x,y,image,xo,yo);

% create noise array with location of NaNs in polar_array
polar_noise = zeros(size(polar_array));
coords = find(isnan(polar_array));
polar_noise(coords) = 1;

polar_array = double(polar_array)./255;

% start diagnostics, writing out eye image with rings overlayed

% get rid of outling points in order to write out the circular pattern
coords = find(xo > size(image,2));
xo(coords) = size(image,2);
coords = find(xo < 1);
xo(coords) = 1;

coords = find(yo > size(image,1));
yo(coords) = size(image,1);
coords = find(yo<1);
yo(coords) = 1;

xo = round(xo);
yo = round(yo);

xo = int32(xo);
yo = int32(yo);

ind1 = sub2ind(size(image),double(yo),double(xo));

image = uint8(image);

image(ind1) = 255;
%get pixel coords for circle around iris
[x,y] = circlecoords([x_iris,y_iris],r_iris,size(image));
ind2 = sub2ind(size(image),double(y),double(x));
%get pixel coords for circle around pupil
[xp,yp] = circlecoords([x_pupil,y_pupil],r_pupil,size(image));
ind1 = sub2ind(size(image),double(yp),double(xp));

image(ind2) = 255;
image(ind1) = 255;

% % write out rings overlaying original iris image
% w = cd;
% % cd(DIAGPATH);
%
% imwrite(image,[eyeimage_filename,'_method',num2str(method),'-normal.jpg'],'jpg');
%
% cd(w);
%
% % end diagnostics

%replace NaNs before performing feature encoding
coords = find(isnan(polar_array));
polar_array2 = polar_array;
polar_array2(coords) = 0.5;
avg = sum(sum(polar_array2)) / (size(polar_array,1)*size(polar_array,2));
polar_array(coords) = avg;

author: Qingbao

repo: iris

save_path: github-repos/MATLAB/Qingbao-iris


text:
%#eml
function R = radix4FFT3_FixPtEML(s)
    % This is a radix-4 FFT, using decimation in frequency
    % The input signal can be floating point or fixed point
    % Works with real or complex input
    % Modified for Embedded MATLAB support

    % Initialize variables and signals
    % NOTE: The length of the input signal should be a power of 4: 4, 16, 64, 256, etc.
    N = length(s);
    M = log2(N)/2;

    % Initialize variables for floating or fixed point sim
    if isfi(s)
        NT = numerictype(s);
        FM = fimath(s);
        wl = NT.WordLength;
        W=fi(exp(-1*2j*pi*(0:N-1)/N), 1, wl, wl - 2, FM);
        S = fi(complex(zeros(size(s))), NT, FM);
        R = fi(complex(zeros(size(s))), 1, wl, wl  - 1 - 2*M, FM);
        sTemp = fi(complex(zeros(size(s))), NT, FM);
        x = fi(complex(zeros(1,4)), NT, FM);
    else
        W=exp(-1*2j*pi*(0:N-1)/N);
        S = complex(zeros(size(s)));
        R = complex(zeros(size(s)));
        sTemp = complex(zeros(size(s)));
        x = complex(zeros(1,4));
    end

    % FFT algorithm
    % Calculate butterflies for first M-1 stages
    sTemp = s;
    for stage = 0:M-2
        for n=1:N/4
            for m=1:4
                x(m) = sTemp(n + (m-1)*N/4);
            end
            S((1:4)+(n-1)*4) = radix4bfly(x, floor((n-1)/(4^stage)) *(4^stage), 1, W);
        end
        sTemp = S;
    end

    % Calculate butterflies for last stage
    for n=1:N/4
        for m=1:4
            x(m) = sTemp(n + (m-1)*N/4);
        end
        S((1:4)+(n-1)*4) = radix4bfly(x, floor((n-1)/(4^stage)) * (4^stage), 0, W);
    end

    % Rescale the final output
    R(:) = S*N;

end

function Z = radix4bfly(x,segment,stageFlag,W)
    % For the last stage of a radix-4 FFT all the ABCD multiplers are 1.
    % Use the stageFlag variable to indicate the last stage
    % stageFlag = 0 indicates last FFT stage, set to 1 otherwise

    % Initialize variables and scale by 1/4
    if isfi(x)
        a=bitshift(x(1),-2);b=bitshift(x(2),-2);c=bitshift(x(3),-2);d=bitshift(x(4),-2);
    else
        a=x(1)*.25;b=x(2)*.25;c=x(3)*.25;d=x(4)*.25;
    end

    % Radix-4 Algorithm
    A=a+b+c+d;
    B=(a-b+c-d)*W(2*segment*stageFlag + 1);
    C=(a-complex(-imag(b),real(b))-c+complex(-imag(d),real(d)))*W(segment*stageFlag + 1);
    D=(a+complex(-imag(b),real(b))-c-complex(-imag(d),real(d)))*W(3*segment*stageFlag + 1);

    Z = [A B C D];

end


text:
function vis_sticks(expidxs,nImgs,bSort,bTrain)

fprintf('vis_sticks()\n');

if (nargin < 3)
    bSort = false;
end

if (nargin < 4)
    bTrain = false;
end

nJoints = 14;

for expidx = expidxs
    fprintf('**********************************************************************************\n');
    % load experiment parameters
    p = exp_params(expidx);
    % load ground truth
    if (bTrain)
        load(p.trainGT);
        path = fileparts(p.evalTrain);
    else
        load(p.testGT);
        path = fileparts(p.evalTest);
    end
    if (~exist('annolist','var'))
        annolist = single_person_annolist;
    end
    annolist = set_full_path(annolist,p.imgDir);
    fnameDist = [path '/distAll'];

    try
%         assert(false);
        load(fnameDist,'pidxs','keypointsAll','matchPCK');
        assert(exist('keypointsAll','var') && length(pidxs) == length(p.pidxs) && ...
               sum(pidxs ~= p.pidxs) == 0 && size(matchPCK,2) == nJoints);
    catch
        warning('keypoints are not found, running evalPCK...');
        evalPCK(expidx);
        load(fnameDist,'pidxs','keypointsAll','matchPCK');
        assert(exist('keypointsAll','var') && length(pidxs) == length(p.pidxs) && ...
               sum(pidxs ~= p.pidxs) == 0 && size(matchPCK,2) == nJoints);
    end

    nImgs = min(length(keypointsAll),nImgs);
    keypointsAll = keypointsAll(1:nImgs);
    if (isfield(p, 'testGTnopad'))
        load(p.testGTnopad,'annolist');
        keypointsAll = projectKeypoints(keypointsAll,fileparts(p.testGTnopad),fileparts(p.testGT));%,annolist
        endpointsAllGT = gt2endpoints(p.testGTnopad);
    else
        endpointsAllGT = gt2endpoints(p.testGT);
    end
    endpointsAll = keypoints2endpoints(keypointsAll);
    matchPCP = eval_match_parts_gt(endpointsAll, annolist(1:min(length(annolist),length(endpointsAll))), 1.0, 1);
    pcp = eval_plot_pcp(matchPCP,1:10);

    if (bSort)
        visDir = [p.expDir '/' p.shortName '/cachedir/test/vis/sort'];
        pckImg = zeros(nImgs,1);
        nanImg = zeros(nImgs,1);
        for imgidx = 1:nImgs
            pckImg(imgidx) = sum(matchPCK(imgidx,~isnan(matchPCK(imgidx,:))));
            allNaNImg(imgidx) = sum(~isnan(matchPCK(imgidx,:))) == 0;
        end
        pckImg(allNaNImg) = nan;
        [~,idxs] = sort(pckImg);
        endpointsAll = endpointsAll(idxs,:);
        keypointsAll = keypointsAll(idxs,:);
        matchPCP = matchPCP(idxs,:);
        matchPCK = matchPCK(idxs,:);
        annolist = annolist(idxs);
        endpointsAllGT = endpointsAllGT(idxs);
    else
        visDir = [p.expDir '/' p.shortName '/cachedir/test/vis/'];
    end
    plot_sticks(p, endpointsAll, annolist, visDir, matchPCP, keypointsAll, matchPCK, p.pidxs, endpointsAllGT, bSort);
end

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [resCe, rhsCe] = electrolyteDiffusion(ce,dCe,jflux,T,param)
% electrolyteDiffusion evaluates the residual for the electrolyte
% concentration of Li-ions in the electrolyte solution.

%   This file is part of the LIONSIMBA Toolbox
%
%	Official web-site: 	http://sisdin.unipv.it/labsisdin/lionsimba.php
% 	Official GitHUB: 	https://github.com/lionsimbatoolbox/LIONSIMBA
%
%   LIONSIMBA: A Matlab framework based on a finite volume model suitable for Li-ion battery design, simulation, and control
%   Copyright (C) 2016-2018 :Marcello Torchio, Lalo Magni, Davide Raimondo,
%                            University of Pavia, 27100, Pavia, Italy
%                            Bhushan Gopaluni, Univ. of British Columbia,
%                            Vancouver, BC V6T 1Z3, Canada
%                            Richard D. Braatz,
%                            Massachusetts Institute of Technology,
%                            Cambridge, Massachusetts 02142, USA
%
%   Main code contributors to LIONSIMBA 2.0:
%                           Ian Campbell, Krishnakumar Gopalakrishnan,
%                           Imperial college London, London, UK
%
%   LIONSIMBA is a free Matlab-based software distributed with an MIT
%   license.

% Diffusion coefficients
% Comment this for benchmark purposes
Deff_p = param.ElectrolyteDiffusionFunction(ce(1:param.Np),T(param.Nal+1:param.Nal+param.Np),param,'p');
Deff_s = param.ElectrolyteDiffusionFunction(ce(param.Np+1:param.Np+param.Ns),T(param.Nal+param.Np+1:param.Nal+param.Np+param.Ns),param,'s');
Deff_n = param.ElectrolyteDiffusionFunction(ce(param.Np+param.Ns+1:end),T(param.Nal+param.Np+param.Ns+1:end-param.Ncu),param,'n');

% Uncomment this for benchmark purposes
% Deff_p = repmat(param.Dp*param.eps_p^param.brugg_p,param.Np,1);
% Deff_s = repmat(param.Ds*param.eps_s^param.brugg_s,param.Ns,1);
% Deff_n = repmat(param.Dn*param.eps_n^param.brugg_n,param.Nn,1);

% Interpolation of the diffusion coefficients
[Deff_p, Deff_s, Deff_n] = interpolateDiffusionCoefficients(Deff_p,Deff_s,Deff_n,param);

%% Positive A Matrix
A_p = - diag(Deff_p);

A_p(2:end,2:end)    = A_p(2:end,2:end) - diag(Deff_p(1:end-1));
A_p(1:end-1,2:end)  = A_p(1:end-1,2:end) + diag(Deff_p(1:end-1));
A_p(2:end,1:end-1)  = A_p(2:end,1:end-1) + diag(Deff_p(1:end-1));
%% Separator A Matrix
A_s = - diag(Deff_s);

A_s(2:end,2:end)    = A_s(2:end,2:end) - diag(Deff_s(1:end-1));
A_s(1:end-1,2:end)  = A_s(1:end-1,2:end) + diag(Deff_s(1:end-1));
A_s(2:end,1:end-1)  = A_s(2:end,1:end-1) + diag(Deff_s(1:end-1));
%% Negative A Matrix
A_n = - diag(Deff_n);

A_n(2:end,2:end)    = A_n(2:end,2:end) - diag(Deff_n(1:end-1));
A_n(1:end-1,2:end)  = A_n(1:end-1,2:end) + diag(Deff_n(1:end-1));
A_n(2:end,1:end-1)  = A_n(2:end,1:end-1) + diag(Deff_n(1:end-1));
% Fix the last elements of the A_n
A_n(end,end-1:end)  = [Deff_n(end-1) -Deff_n(end-1)];
%% A_tot matrix
A_tot = blockDiagonalMatrix(param,A_p,A_s,A_n);

% Divide by the deltax and the length of the positive electrode
A_tot(1:param.Np,1:param.Np)    = A_tot(1:param.Np,1:param.Np)/(param.deltax_p^2*param.len_p^2);
% Reset values on the lines for the interfaces conditions
A_tot(param.Np,:)           = 0;
A_tot(param.Np+1,:)         = 0;

% Divide by the deltax and the length of the separator
A_tot(param.Np+1:param.Np+param.Ns,param.Np+1:param.Np+param.Ns)    = A_tot(param.Np+1:param.Np+param.Ns,param.Np+1:param.Np+param.Ns)/(param.deltax_s^2 * param.len_s^2);
% Reset values on the lines for the interfaces conditions
A_tot(param.Np+param.Ns,:)      = 0;
A_tot(param.Np+param.Ns+1,:)    = 0;

% Divide by the deltax and the length of the negative electrode
A_tot(param.Np+param.Ns+1:end,param.Np+param.Ns+1:end)      = A_tot(param.Np+param.Ns+1:end,param.Np+param.Ns+1:end)/(param.deltax_n^2 * param.len_n^2);

%% Interface between separator and positive electrode (last volume in the positive electrode)

% Compute the common denominator at the interface
den_s   = (param.deltax_p*param.len_p/2 + param.deltax_s*param.len_s/2);
% Last diffusion coefficient of the positive electrode
last_p  = Deff_p(end-1)/(param.deltax_p*param.len_p);
% Diffusion coefficient on the interface
first_s = Deff_p(end)/den_s;
% Fix the values at the boundaries
A_tot(param.Np,param.Np-1:param.Np+1) = [last_p -(last_p+ first_s) first_s]/(param.deltax_p*param.len_p*param.eps_p);

%% Interface between separator and positive electrode (first volume in the separator)

% Compute the common denominator at the interface
den_s       = (param.deltax_p*param.len_p/2 + param.deltax_s*param.len_s/2);
% First diffusion coefficient in the separator
second_s    = Deff_s(1)/(param.deltax_s*param.len_s);
% Diffusion coefficient on the interface
first_s     = Deff_p(end)/den_s;

A_tot(param.Np+1,param.Np:param.Np+2) = [first_s -(first_s+second_s) second_s]/(param.deltax_s*param.len_s*param.eps_s);

%% Interface between separator and negative electrode (last volume in the separator)

% Compute the common denominator at the interface
den_s   = (param.deltax_s*param.len_s/2 + param.deltax_n*param.len_n/2);
% Last diffusion coefficient in the separator
last_s  = Deff_s(end-1)/(param.deltax_s*param.len_s);
% Diffusion coefficient on the interface
first_n = Deff_s(end)/den_s;

A_tot(param.Np+param.Ns,param.Np+param.Ns-1:param.Np+param.Ns+1) = [last_s -(last_s+first_n) first_n]/(param.deltax_s*param.len_s*param.eps_s);

%% Interface between separator and negative electrode (first volume in the negative electrode)

% Compute the common denominator at the interface
den_n       = (param.deltax_s*param.len_s/2 + param.deltax_n*param.len_n/2);
% First diffusion coefficient in the negative electrode
second_n    = Deff_n(1)/(param.deltax_n*param.len_n);
% Diffusion coefficient on the interface
first_n     = Deff_s(end)/den_n;

A_tot(param.Np+param.Ns+1,param.Np+param.Ns:param.Np+param.Ns+2) = [first_n -(first_n+second_n) second_n]/(param.deltax_n*param.len_n*param.eps_n);

%% Useful stuff
a_tot       = [
    repmat(param.a_i(1),param.Np,1);...
    zeros(param.Ns,1);...
    repmat(param.a_i(3),param.Nn,1)...
    ];

jflux_tot   = [
    jflux(1:param.Np);...
    zeros(param.Ns,1);...
    jflux(param.Np+1:end)...
    ];

eps_tot     = [
    repmat(param.eps_p,param.Np,1);...
    repmat(param.eps_s,param.Ns,1);...
    repmat(param.eps_n,param.Nn,1)
    ];

K = 1./eps_tot;
A_eps = diag(K);

% Build porosities matrix
if(~isa(eps_tot,'casadi.MX') && ~isa(eps_tot,'casadi.SX'))
    A_eps = A_eps + diag(K(1:end-1),1);
    A_eps = A_eps + diag(K(1:end-1),-1);
    A_eps = sparse(A_eps);
else
    A_u_l                   = diag(K(1:end-1));
    A_eps(1:end-1,2:end)    = A_eps(1:end-1,2:end)+A_u_l;
    A_eps(2:end,1:end-1)    = A_eps(2:end,1:end-1)+A_u_l;
end

A_eps(param.Np,param.Np-1:param.Np+1) = 1;
A_eps(param.Np+1,param.Np:param.Np+2) = 1;

A_eps(param.Np+param.Ns,param.Np+param.Ns-1:param.Np+param.Ns+1) = 1;
A_eps(param.Np+param.Ns+1,param.Np+param.Ns:param.Np+param.Ns+2) = 1;

G = A_eps.*A_tot;

% Write the RHS of the equation
rhsCe = (G*ce + K.*(1-param.tplus).*a_tot.*jflux_tot);

% Write the residual of the equation
resCe = dCe - rhsCe;
end

author: lionsimbatoolbox

repo: LIONSIMBA


text:
 function ob = ctranspose(ob)
%function ob = ctranspose(ob)
%	"ctranspose" method for Gtomo2 class

ob.is_transpose = ~ob.is_transpose;

[ob.index1 ob.index2] = deal(ob.index2, ob.index1);
ob.dims = fliplr(ob.dims);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function output = perceptronTrain(input)
%perceptronTrain Training of two-class perceptron.
%   OUTPUT = perceptronTrain(INPUT) training of two-class perceptron. The
%   input and output parameters of are structures defined as follows.
%
%   Parameter input is a structure with the following fields:
%
%   input.X             A matrix of size n-by-np whose columns are
%                       n-dimensional input pattern vectors, and np is
%                       the number of such vectors. For example, for a
%                       set of 100, 3-element patterns, matrix inputs.X
%                       would be of size 3-by-100. The function augments
%                       the patterns by appending a 1 at the end of all
%                       pattern vectors.
%
%   input.r             Class membership vector of size 1-by-np whose
%                       kth value is 1 if the pattern vector in the kth
%                       column of X belongs to class c1; otherwise the
%                       kth value of input.r is -1.
%
%   input.Alpha         Learning rate positive scalar that defaults to
%                       0.5.
%
%   input.W0            Initial augmented weight vector of size
%                       (n + 1)-by-1. Its components default to uniform
%                       random numbers in the range [0,1].
%
%   input.NumEpochs     The number of training epochs. Defaults to 100.
%
%
%   Parameter output is a structure with the following fields.
%
%   output.W            (n + 1)-by-1 weight vector at completion of
%                       training. This may or not be a solution that
%                       separates the two classes, depending on the
%                       number of epochs specified and whether the
%                       classes are linearly separable.
%
%   output.Convergence  True (1) if training converged and false (0)
%                       otherwise.
%
%   output.ActualEpochs The actual number of epochs of training
%                       executed.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% AUGMENT THE PATTERN VECTORS.
% Number of patterns.
np = size(input.X,2);
% Augment the pattern vectors.
X = cat(1,input.X,ones(1,np));

% SET DEFAULTS.
if ~isfield(input,'Alpha')
   input.Alpha = 0.5;
end
if ~isfield(input,'NumEpochs')
   input.NumEpochs = 100;
end
if ~isfield(input,'W0')
   rng('shuffle');
   input.W0 = rand(size(X,1),1);
end

% INITIALIZATION.
% Make sure that input.W0 is a column vector.
input.W0 = input.W0(:);
% Initial values for loop.
w_last = input.W0;
w = input.W0;
output.ActualEpochs = 0;

% ITERATE.
for I = 1:input.NumEpochs
   output.ActualEpochs = output.ActualEpochs + 1;
   output.Convergence = true;
   for J = 1:np
      % An error has been committed for the Jth pattern if the sign of
      % w'*input.X(:,J) does not match the sign of input.r
      if sign(w'*X(:,J)) ~= input.r(J)
         w = w_last + input.r(J)*input.Alpha*X(:,J);
         output.Convergence = false;
         w_last = w;
      end
   end
   % Exit if convergence was not set to false for an entire epoch.
    if output.Convergence
      output.W = w_last;
      output.ActualEpochs = I;
      break
    end
end
output.ActualEpochs = I;
output.W = w_last;

author: dipum

repo: dipum-toolbox

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


text:
function [ y, j, f ] = jed_to_yjf_roman ( jed )

%*****************************************************************************80
%
%% JED_TO_YJF_ROMAN converts a JED to a Roman YJF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 June 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real JED, the Julian Ephemeris Date.
%
%    Output, integer Y, J, real F, the YJF date.
%
  [ y1, m1, d1, f1 ] = jed_to_ymdf_roman ( jed );

  [ y, j, f ] = ymdf_to_yjf_roman ( y1, m1, d1, f1 );

  return
end


text:
function [Ex f u v Alfa EtaInf itrs] = FalknerSkanSolver(b0, b, eps1, eps2, Alfa0, EtaInf0, hObject, trace)
    global SOL;
    global ETA_INF;
    global BETA0;
    global BETA;
    global EX;
    t1 = cputime;
    handles = guidata(hObject);
    Ex_span = 0:0.05:1;
    %
    % Allocation memory
	%
	Alfa = [];
	EtaInf = [];
	J = zeros(2);
	F = zeros(2,1);
	% Input values
	Alfa(1) = Alfa0;
	EtaInf(1) = EtaInf0;
    BETA0 = b0;
    BETA = b;
	% Initial Conditions
	p = eps1*10;
	q = eps2*10;
	k = 0;
    if (trace == 1)
        [msg err] = sprintf('  i      Alfa           Eta Inf         q(alfa, eta Inf)');
        set(handles.Info_LstBx, 'String', msg);
        [msg1 err] = sprintf('|-------------------------------------------------------------');
        msg = strcat(msg, msg1);
        set(handles.Info_LstBx, 'String', msg);
        [msg1 err] = sprintf('|%3d  %1.11f      %1.5f    ', k+1, ...
                         Alfa(k+1),EtaInf(k+1));
        msg = strcat(msg, msg1);
        set(handles.Info_LstBx, 'String', msg);
    end
	while ((abs(p) > eps1) || (abs(q) > eps2))
    	k = k + 1;
    	p = eps1*10;
    	ETA_INF = EtaInf(k);
    	while (abs(p) > eps1)
        	IC = [0 0 Alfa(k)]';
       		[EX SOL] = myode45(@FalknerSkanSys, Ex_span, IC,[],0);
            f = SOL(:,1)';
            u = SOL(:,2)';
            v = SOL(:,3)';
            [nodes c] = size(EX);
            p = u(nodes)-1;
            q = v(nodes);
            IC = [0 0 1]';
            [Ex_ d_dalfa] = myode45(@JacobianAlfa, Ex_span, IC,[],0);
            [nodes1 c] = size(Ex_);
            dpdalfa =  d_dalfa(nodes1,2);
            DeltaAlfa = -p/dpdalfa;
            Alfa(k) = Alfa(k) + DeltaAlfa;
        end
        IC = [0 0 0]';
        [Ex_ d_dEtaInf] = myode45(@JacobianEtaInf, Ex_span,IC,[],0);
        [nodes2 c] = size(Ex_);
        % Jacobian matrix
        J (1,1) = d_dalfa(nodes1, 2);
        J (1,2) = d_dEtaInf(nodes2, 2);
        J (2,1) = d_dalfa(nodes1, 3);
        J (2,2) = d_dEtaInf(nodes2, 3);
        F(1) = p;
        F(2) = q;
        Delta = -J\F;
        Alfa(k+1) = Alfa(k)+Delta(1);
        EtaInf(k+1) = EtaInf(k)+Delta(2);
        if (trace == 1)
             [msg1 err] = sprintf('|%3d  %1.11f      %1.5f          %1.4e', k+1, Alfa(k+1),EtaInf(k+1), abs(q) );
             msg = strcat(msg, msg1);
             set(handles.Info_LstBx, 'String', msg);
        end
    end
    % output values
    Ex = EX;
    itrs = k+1;
    t2 = cputime - t1;
    if (trace == 1)
        [msg1 err] = sprintf('| ');
        msg = strcat(msg, msg1);
        set(handles.Info_LstBx, 'String', msg);
        [msg1 err] = sprintf('|Computational time: %d s', t2);
        msg = strcat(msg, msg1);
        set(handles.Info_LstBx, 'String', msg);
    end


text:
function t = uminus(t)
%UMINUS Unary minus (-) for tensors.
%
%   See also TENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

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

t.data = -t.data;

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt