text:
 function xs = eml_osps(x, Gb, yi, ci, ri, niter, pixmax, curv, relax0, chat)
%function xs = eml_osps(x, Gb, yi, ci, ri, niter, pixmax, curv, relax0, chat)
% E-ML-OSPS algorithm for emission Poisson problem
% (ordered subsets separable paraboloidal surrogates)
% model: Y_i ~ Poisson(c_i [G x]_i + r_i)
% in
%	x	[np,1]		initial estimate
%	Gb			Gblock object (see eml_osps_test.m)
%	yi,ci,ri [nb,na]	see em_fbp.m (for model too)
%	niter			# iterations
%	pixmax			upper constraint for pixel values
%	curv	'oc' for erdogan's optimal curvatures
%		'pc' for erdogan's fast precomputed curvatures,
%			which usually provides faster convergence,
%			but can be nonmonotone
%	relax0	[1] or [2]	relax0 or (relax0, relax_rate)
%	chat
% out
%	xs [np,niter]	updated image vectors each iteration
%
% Copyright Mar 2000, Jeff Fessler, The University of Michigan

if nargin < 3, ir_usage, end

nblock = block_ob(Gb, 'n');
starts = subset_start(nblock);

if ~isvar('ci') || isempty(ci)
	ci = ones(size(yi));
end
if ~isvar('ri') || isempty(ri)
	ri = zeros(size(yi));
end
if ~isvar('niter') || isempty(niter),	niter = 1;	end
if ~isvar('pixmax') || isempty(pixmax),	pixmax = inf;	end
if ~isvar('curv') || isempty(curv),	curv = 'oc';	end
if ~isvar('chat') || isempty(chat),	chat = true;	end

if ~isvar('relax0') || isempty(relax0)
	relax0 = 1;
end
if length(relax0) == 1
	relax_rate = 0;
elseif length(relax0) == 2
	relax_rate = relax0(2);
	relax0 = relax0(1);
else
	error relax
end

eml_check(yi, ci, ri);

[nb, na] = size(yi);

gi = sum(Gb')';		% g_i = sum_j g_ij
gi = reshape(gi, nb, na);

%
%	precomputed curvatures
%
denom = zeros(numel(x), nblock);
if streq(curv, 'pc')
	ni = eml_curvature(yi, ci, ri, [], [], curv);
	denom = Gb' * col(gi .* ni);
%	printf('ni range %g %g', min(ni(:)), max(ni(:)))
%	printf('denom range %g %g', min(denom(:)), max(denom(:)))
end

%
%	loop over iterations
%
xs = zeros(numel(x), niter);
x = max(x,0);
x = min(x,pixmax);
xs(:,1) = x;
for iter = 2:niter

	relax = 1;

	%
	% loop over subsets
	%
	for iset=1:nblock
		iblock = starts(iset);
		ia = iblock:nblock:na;

		li = Gb{iblock} * x;			% l=G*x "line integrals"
		li = reshape(li, nb, length(ia));
		yb = ci(:,ia) .* li + ri(:,ia);		% predicted meas. means

		% fix: need to be careful here with 0/0 -> 0
		dothi = ci(:,ia) .* (yi(:,ia) ./ yb - 1);

		% non-precomputed curvatures (notably, optimal curvature),
		% for ensured monotone increase
		if ~streq(curv, 'pc')
			ni = eml_curvature(yi(:,ia), ci(:,ia), ri(:,ia), li, yb, curv);
			denom = nblock * (Gb{iblock}' * col(gi(:,ia) .* ni));
		end

		grad = Gb{iblock}' * dothi(:);
		num = nblock * grad;
		x = x + relax * num ./ denom;	% relaxed update
		x = max(x,0);			% lower bound
		x = min(x,pixmax);		% upper bound
	end

	xs(:,iter) = x;
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [h]=rayleigh(fd,t)

    N=40;
    wm=2*pi*fd;
    N0=N/4;
    Tc=zeros(1,length(t));
    Ts=zeros(1,length(t));
    P_nor=sqrt(1/N0);
    theta=2*pi*rand(1,1)-pi;

    for ii=1:N0

            alfa(ii)=(2*pi*ii-pi+theta)/N;
            fi_tc=2*pi*rand(1,1)-pi;
            fi_ts=2*pi*rand(1,1)-pi;
            Tc=Tc+cos(cos(alfa(ii))*wm*t+fi_tc);
            Ts=Ts+cos(sin(alfa(ii))*wm*t+fi_ts);
    end;

   h=P_nor*(Tc+j*Ts );

author: liyanluminary


text:
function [safety, order, union, alg_conn, safety_obs, min_d_obs ] = ...
    compute_swarm_performance(pos_history, vel_history, ...
    p_swarm, dirname)
% Compute swarm performance - This function allows to compute the
% performance from the history of the swarming variables (time, position,
% velocity, acceleration).
%
% Inputs:
%   pos_history: series of agents' positions
%   vel_history: series of agents' velocities
%   p_swarm: structure with swarm parameters
%   dirname: path of the directory where results are saved
%
% Outputs:
%   safety: safety metric. It reflects the number of agent-agent collisions
%   order: order metric. It measures the correlation between velocity
%               vectors
%   union: union metric. It measures the number of connected components of
%               the undirected graph
%   alg_conn: algebraic connectivity metric.
%   safety_obs: safety againt obstacles metric. reflects the number of
%               agent-obstacle collisions
%   min_d_obs: minimum distance agent-obstacle
%

%% Init variables

[t_steps,nx] = size(pos_history);
nb_agents = nx/3;
M = zeros(t_steps,nb_agents,nb_agents);
nb_ag_coll = zeros(t_steps,1);
nb_obs_coll = zeros(t_steps,1);
nb_conn_comp = zeros(t_steps,1);

safety = ones(t_steps,1);
order = zeros(t_steps,1);
union = zeros(t_steps,1);
alg_conn = zeros(t_steps,1);
safety_obs = zeros(t_steps,1);
min_d_obs = zeros(t_steps,nb_agents);

%% Loop over time

for k = 1:t_steps

    %% Safety: reflects the number of collisions among the swarm agents

    pos_k = pos_history(k,:);
    pos_k = reshape(pos_k,3,[]);
    dist_vect_k = pdist(pos_k');
    if ~isempty(dist_vect_k)
        nb_ag_coll(k) = sum(dist_vect_k < 2*p_swarm.r_coll);
        safety(k) = 1 - (sum(nb_ag_coll(k)) / length(dist_vect_k));
    else
        nb_ag_coll(k)=0;
        safety(k) = 1;
    end

    %% Order: reflects the correlation of the velocity vectors
    distance_matrix_k = squareform(dist_vect_k);
    M(k,:,:) =  compute_neighborhood(distance_matrix_k, p_swarm.r, p_swarm.max_neig);
    for agent = 1:nb_agents
        neig = (M(k,:,agent)==true);
        nn = sum(neig);
        vel_k = vel_history(k,:);
        vel_k = reshape(vel_k,3,[]);
        if nn ~= 0
            scalar_prod = vel_k(:, agent)' * vel_k(:, neig);
            speed_agent = norm(vel_k(:, agent));
            speed_neig = sqrt(sum((vel_k(:, neig) .^ 2), 1));
            order(k) = order(k) + ...
                sum (scalar_prod ./ (speed_agent * speed_neig)) / nn;
        else
            order(k) = order(k) + 1;
        end
    end
    order(k) = order(k)/nb_agents;

    %% Union: reflects the number of connected components in the related
    % symeetric graph

    M_k = squeeze(M(k,:,:));
    A = ((M_k + M_k') > 0);
    [nb_conn_comp(k), ~, ~]  = network_components(A);
    union(k) = (nb_agents - nb_conn_comp(k))/(nb_agents-1);

    %% Connectivity: reflects the algebraic connectivity in the related
    % symmetric graph

    alg_conn(k) = compute_alg_connectivity(A)/nb_agents;

    %% Safety with obstacles: reflects the number of collisions between swarm agents and obstacles

    nb_possible_coll = 0;
    if p_swarm.is_active_spheres | p_swarm.is_active_cyl | p_swarm.is_active_arena

        pos_k = pos_history(k,:);
        pos_k = reshape(pos_k,3,[]);

        if p_swarm.is_active_spheres
            c_spheres = p_swarm.spheres(1:3,:);
            r_spheres = p_swarm.spheres(4,:);

            D_spheres = pdist2(pos_k',c_spheres');
            min_d_obs(k,:) = min(pdist2(pos_k(1:2,:)',c_spheres') - repmat(r_spheres, nb_agents, 1),[],2);
            nb_obs_coll(k) = nb_obs_coll(k) + sum(sum(D_spheres < repmat(r_spheres, nb_agents, 1)));
            nb_possible_coll = nb_possible_coll + nb_agents*length(r_spheres);
        end
        if p_swarm.is_active_cyl
            c_cyl = p_swarm.cylinders(1:2,:);
            r_cyl = p_swarm.cylinders(3,:);

            D_cyl = pdist2(pos_k(1:2,:)',c_cyl');
            min_d_obs(k,:) = min(pdist2(pos_k(1:2,:)',c_cyl') - repmat(r_cyl, nb_agents, 1),[],2);
            nb_obs_coll(k) = nb_obs_coll(k) + sum(sum(D_cyl < repmat(r_cyl, nb_agents, 1)));
            nb_possible_coll = nb_possible_coll + nb_agents*length(r_cyl);
        end
        if p_swarm.is_active_arena
            x_coll = sum(pos_k(1,:)< p_swarm.x_arena(1,1) | pos_k(1,:)> p_swarm.x_arena(1,2));
            y_coll = sum(pos_k(2,:)< p_swarm.x_arena(2,1) | pos_k(2,:)> p_swarm.x_arena(2,2));
            z_coll = sum(pos_k(3,:)< p_swarm.x_arena(3,1) | pos_k(3,:)> p_swarm.x_arena(3,2));
            nb_obs_coll(k) = nb_obs_coll(k) + x_coll + y_coll + z_coll;
            nb_possible_coll = nb_possible_coll + nb_agents;
        end

        safety_obs(k) = 1 - (nb_obs_coll(k)/nb_possible_coll);
    end

end

%% Save workspace

if ~isempty(dirname)
    path = strcat(dirname,'/performance');
    save(path);
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function connectivity = compute_alg_connectivity(A)
% Compute alg connectivity - compute algebraic connectivity of the graph
% associated to the swarm at a given time step.

D = diag(sum(A, 2)); % degree matrix
L = D - A;           % laplacian matrix
eigenvalues = eig(L);
eigenvalues = sort(eigenvalues);
connectivity = eigenvalues(2);

end

author: lis-epfl

repo: swarmlab

save_path: github-repos/MATLAB/lis-epfl-swarmlab


text:
function [mx,my,mz] = blochsim2(Mi, bx, by, bz, T1, T2, dt)

%	Evolve magnetization field using time-discretized Bloch equation
%	Input
%		Mi [3,dim]	initial X,Y,Z magnetization
%					Note: normalized to magnitude <= Mo=1
%		bx,by,bz [ntime,dim]	effective X,Y,Z applied magnetic field
%					(Tesla), for rotating frame (no Bo)
%		T1	[dim]		spin-lattice relaxation time (msec)
%		T2	[dim]		spin-spin relaxation time (msec)
%		dt	scalar		time interval (msec)
%	Output
%		mx,my,mz [ntime,dim]  X,Y,Z magnetization as a function of time
%
% uses:
% https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula

% constants
gambar = 42.57e3;        % gamma/2pi in kHz/T
gam = 267513; %gambar*2*pi;
% Put Beff into units rotations, T1 and T2 into losses
bx = bx.*(dt*gam);      % rotation angle/step
by = by.*(dt*gam);      % rotation angle/step
bz = bz.*(dt*gam);      % rotation angle/step
% Put relaxations into losses/recovery per step
T1 = (dt ./ T1(:));
T2 = (1 - dt ./ T2(:));

% size checks
if ~isreal(bx) || ~isreal(by) || ~isreal(bz)
    bx = real(bx);
    by = real(by);
    bz = real(bz);
    disp('Warning: B field must be real valued - using only the real part');
end
if (size(bx) ~= size(by)) || (size(bx) ~= size(bz))
    disp('Error: B vectors not the same length')
    return;
end
if (size(Mi,2) ~= size(bx,2))
    disp('Error: Initial magnetization not right size')
    return;
end
if (size(Mi,2) ~= length(T1))
    disp('Error: T1 vector not right size')
    return;
end
if (size(Mi,2) ~= length(T2))
    disp('Error: T2 vector not right size')
    return;
end

nstep = size(bx,1);
%
% Initialize outputs
mx = zeros(size(bx,1)+1,size(bx,2)); % record one more location
my = zeros(size(mx));
mz = zeros(size(mx));
mx(1,:) = Mi(1,:);
my(1,:) = Mi(2,:);
mz(1,:) = Mi(3,:);

% stable bloch equation simulator: rotations are explicitly
% calculated and carried out on the magnetization vector
for lp = 2:nstep+1 % Hao: modified from 2:nstep to 2:nstep+1. This will record the mag. after last pulse point.
  B = [bx(lp-1,:); by(lp-1,:); bz(lp-1,:)]';
  %	Compute sines & cosines of field angles:
  %	Theta = angle w.r.t positive z axis
  %	Phi   = angle w.r.t positive x axis
  %	Psi   = angle w.r.t transformed positive x axis
  %
  Bmag = sqrt(sum(B.^2,2));		% Magnitude of applied field
  Btrans = sqrt(B(:,1).^2 + B(:,2).^2);	% Magnitude of transverse applied field
  ct = ones(size(B,1),1);
  good = Bmag ~= 0;
  if any(good)
	ct(good) = B(good,3) ./ Bmag(good);	% cos(theta)
  end
  st = sqrt(1 - ct.^2);				% sin(theta) > 0

  cphi = ones(size(B,1),1);
  good = Btrans ~= 0;
  if any(good)
	cphi(good) = B(good,1) ./ Btrans(good);	% cos(phi)
  end
  sphi = sqrt(1 - cphi.^2) .* sign(B(:,2));	% sin(phi)

  cpsi = cos(Bmag);			% cos(psi)
  spsi = sin(Bmag);			% sin(psi)

  %
  %	Evolve
  %
  if any(Bmag ~= 0)
    Mx0 = mx(lp-1,:)';
    My0 = my(lp-1,:)';
    Mz0 = mz(lp-1,:)';

    Mx1 = cphi.*(ct.*(cpsi.*(ct.*(sphi.*My0+cphi.*Mx0)-st.*Mz0) ...
    + spsi.*(cphi.*My0-sphi.*Mx0))+st.*(ct.*Mz0+st.*(sphi.*My0+cphi.*Mx0))) ...
    - sphi.*(-spsi.*(ct.*(sphi.*My0+cphi.*Mx0)-st.*Mz0) ...
    + cpsi.*(cphi.*My0-sphi.*Mx0));
    My1 = sphi.*(ct.*(cpsi.*(ct.*(sphi.*My0+cphi.*Mx0)-st.*Mz0) ...
    + spsi.*(cphi.*My0-sphi.*Mx0))+st.*(ct.*Mz0+st.*(sphi.*My0+cphi.*Mx0))) ...
    + cphi.*(-spsi.*(ct.*(sphi.*My0+cphi.*Mx0)-st.*Mz0) ...
    + cpsi.*(cphi.*My0-sphi.*Mx0));
    Mz1 = ct.*(ct.*Mz0+st.*(sphi.*My0+cphi.*Mx0)) ...
    - st.*(cpsi.*(ct.*(sphi.*My0+cphi.*Mx0)-st.*Mz0) ...
    + spsi.*(cphi.*My0-sphi.*Mx0));
  else
    Mx1 = mx(lp-1,:)';
    My1 = my(lp-1,:)';
    Mz1 = mz(lp-1,:)';
  end
  % relaxation effects: "1" in Mz since Mo=1 by assumption
  mx(lp,:) = (Mx1 .* T2)';
  my(lp,:) = (My1 .* T2)';
  mz(lp,:) = (Mz1 + (1- Mz1).* T1)'; % true, since mz = Mo - (Mo-mz)*exp(-dt/T1) ~ Mo - (Mo-mz)(1-dt/T1) = mz + (Mo-mz)*dt/T1

end % end loop through time

% Hao: made the change so my,my,mz will NOT record the initial mag. But record the mag. after last pulse point
  mx(1,:) = [];
  my(1,:) = [];
  mz(1,:) = [];   % true, since mz = Mo - (Mo-mz)*exp(-dt/T1) ~ Mo - (Mo-mz)(1-dt/T1) = mz + (Mo-mz)*dt/T1


text:
clear
clc
A=[2 1 5 3 4];
B=[1 5 3 4 2];
C=[4 1 2 3 5];
D=[1 4 2 3 5];
E=[5 4 3 1 2];
neighbork=[A;B;C;D;E];

center_route=Center(neighbork);

author: sfvsfv

repo: Mathematical-modeling


text:
function l = s_to_l ( s )

%*****************************************************************************80
%
%% S_TO_L reads a logical value from a string.
%
%  Discussion:
%
%    There are several ways of representing logical data that this routine
%    recognizes:
%
%      False   True
%      -----   ----
%
%      0       1
%      F       T
%      f       t
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be read.
%
%    Output, integer L, the logical value read from the string.
%
  s_length = length ( s );

  for i = 1 : s_length
    if ( s(i) == '0' || s(i) == 'F' || s(i) == 'f' )
      l = 0;
      return
    elseif ( s(i) == '1' || s(i) == 'T' || s(i) == 't' )
      l = 1;
      return
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'S_TO_L - Fatal error!\n' );
  fprintf ( 1, '  The input string did not contain logical data.\n' );
  error ( 'S_TO_L - Fatal error!' );
end

author: johannesgerer

repo: jburkardt-m


text:
classdef MaterialFactory < handle

    methods (Access = public, Static)

        function material = create(cParams)
            switch cParams.ptype
                case {'ELASTIC','DIFF-REACT'}
                    switch cParams.pdim
                        case '2D'
                            cParams.nstre = 3;
                            material = Isotropic2dElasticMaterial(cParams);
                        case '3D'
                            cParams.nstre = 6;
                            material = Isotropic3dElasticMaterial(cParams);
                    end

                case 'HYPERELASTIC'
                    switch cParams.pdim
                        case '2D'
                            cParams.connec = cParams.mesh.connec;
                            cParams.dNdx   = cParams.geometry.dNdx;
                            cParams.nnode  = cParams.mesh.nnodeElem;
                            cParams.coord  = cParams.mesh.coord;
                            material = Isotropic2dHyperElasticMaterial(cParams);
                    end

                case 'THERMAL'
                    error('Still not implemented.')
                case 'Stokes'
                    material = Material_Stokes(cParams);
                otherwise
                    error('Invalid ptype.')
            end

        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function el = export_referenced_structure_set_sequence(args)
%"export_referenced_structure_set_sequence"
%   Subfunction to handle referenced_structure_set sequences within the
%   rt_dvh module.  Uses the same layout and principle as the
%   parent function.
%
%   This function takes a CERR DVHs element.
%
%JRA 07/10/06
%NAV 07/19/16 updated to dcm4che3
%   replaced ml2dcm_Element to data2dcmElement
%
%Usage:
%   @export_referenced_structure_set_sequence(args)
%
% 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/>.

%Init output element to empty.
el = [];

%Unpack input data.
tag         = args.tag;
DVHS        = args.data{1};
template    = args.template;

switch tag
    case 528720     %0008,1150  Referneced SOP Class UID
        data = DVHS.Referenced_Structure_Set_SOP_Class_UID;
        el = data2dcmElement(data, tag);

    case 528725     %0008,1155  Referenced SOP Instance UID
        data = DVHS.Referenced_Structure_Set_SOP_Instance_UID;
        el = data2dcmElement(data, tag);

    otherwise
        warning(['No methods exist to populate DICOM ROI_contour module''s ROI_contour_sequence field: ' dec2hex(tag,8) '.']);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [dateNum, dateStr] = gitLastRevDate(repo, file)
dateNum = 0;
dateStr = '';

if ~exist('repo','var') || isempty(repo)
    repo = pwd;
end
if ~exist('file','var') || isempty(file)
    file = '';
end
repoFull = filesepStandard_startup(repo, 'full');
currdir = pwd;

if ismac()
    cmdprefix = 'TERM=ansi; ';
    e = system(cmdprefix);
    if e ~= 0
        return;
    end
else
    cmdprefix = '';
end

ii = 1;
cmds{ii,1} = sprintf('cd %s', repoFull); ii = ii+1;
cmds{ii,1} = sprintf('%sgit log -1 --format=%%ci --date=short %s', cmdprefix, file);
[errs, msgs] = exeShellCmds(cmds, false, true);
if all(errs==0)
    if ~isempty(msgs{2})
        ds = msgs{2};
    else
        ds = '';
        %ds = char(datetime('now','TimeZone','local','Format','yyyy-MM-dd HH:mm:ss'));
    end
    [dateNum, dateStr] = datestr2datenum(ds);
end
cd(currdir)

author: BUNPC


text:
function twod_to_vtk_test ( )

%*****************************************************************************80
%
%% TWOD_TO_VTK_TEST tests the TWOD_TO_VTK library.
%
%  Discussion:
%
%    The data represents a flow solution in a 10 by 3 channel.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 December 2010
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( '\n' );
  fprintf ( 1, 'TWOD_TO_VTK_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the TWOD_TO_VTK library.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read in TWOD finite element data.\n' );
  fprintf ( 1, '  Write it out as a legacy VTK file.\n' );
%
%  Load the data from text files.
%
  nodes = load ( 'channel_nodes.txt' );
  elements = load ( 'channel_elements.txt' );
  uvp = load ( 'channel_values.txt' );

  output_filename = 'channel.vtk';

  title = 'UVP for 10x3 2D channel problem';
%
%  Have the data written to a file.
%
  twod_to_vtk ( nodes, elements, uvp, output_filename, title );
%
%  Terminate.
%
  fprintf ( '\n' );
  fprintf ( 1, 'TWOD_TO_VTK_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( '\n' );
  timestamp ( );

  return
end


text:
function sR=icassoProjection(sR, method, varargin)
%function sR=icassoProjection(sR,[method],['identifier1',val1,'identifier2',val2,...]))
%
%PURPOSE
%
%To project points on plane so that Euclidean distances between the
%projected points correspond to the similarity matrix between IC
%estimates in the Icasso result structure.
%
%EXAMPLES OF BASIC USAGE
%
%   sR=icassoProjection(sR);
%
%makes a CCA projection using default parameters. This is equivalent
%of giving command:
%
%   sR=icassoProjection(sR,'cca', 's2d','sim2dis2','epochs',70,'alpha',0.7);
%
%In the following example 'sim2dis' (i.e., D=1-S) is used to make
%the similarity-to-dissimilarity transformation and a longer,
%sammon's projection is used and a longer training sequence is engaged:
%
%   sR=icassoProjection(sR,'sammon','s2d','sim2dis','epochs',200):
%
%INPUT
%
%[An argument in brackets is optional. If it isn't  given or it's
% an empty matrix/string, the function will use a default value.]
%
% sR       (struct) Icasso result data structure
% [method] (string) 'cca' (defalut) | 'mmds' | 'sammon'
%
%Optional input arguments are given as argument identifier - value
%pairs: 'identifier1', value1, 'identifier2', value2,...
%(case insensitive)
%
%  'epochs'     (scalar) see details
%  'alpha',     (scalar) see details
%  'radius'     (scalar) see details
%  's2d'        (string) see details
%
%OUTPUT
%
% sR (struct) updated Icasso result data structure,
%
%The function updates _only_ the following fields:
%sR.projection.method, sR.projection.parameters, and
%sR.projection.coordinates. See icassoStruct.
%
%DETAILS
%
%The function transforms the similarities S in field
%sR.cluster.similarities into dissimilarities using function
%sim2dis2 as default. Note that this may be different from
%the transformation that was used for making clustering.
%Input argument pair 's2d',<string> gives the name of function
%that is used to make the transformation from similarities
%S=sR.cluster.similarity. There are two ready made functions
%sim2dis.m and sim2dis2.m that can be used:
%   function      makes transformation
%  'sim2dis2'     D=sqrt(1-S)   (default)
%  'sim2dis'      D=1-S
%or you can specify your own function.
%
%The function can use three methods to do the projection on D
%1. Curvilinear Component Analysis (CCA) (preferred)
%2. Principal Coordinates (linear Metric Multi-Dimensional Scaling,MMDS)
%3. Sammon's projection (Sammon)
%CCA and Sammon require some parameters. The default values are set
%according to experience and can be altered if necessary. MMDS is
%automatically used as an initial projection for both Sammon and CCA.
%
%'epochs'  (scalar) the number of epochs for CCA or Sammon.
%           Ignored for MMDS.
%           default(s): CCA: 75, Sammon: 100
%
%'radius'  (scalar) CCA initial radius. Ignored for Sammon and MMDS
%            default(s) CCA: max(10,M/20) where M is the number of
%            estimates
%
%'alpha'   (scalar) learning rate factor for CCA or Sammon.
%            ignored for MMDS
%            default(s): CCA: 0.7, Sammon: 0.7
%
%SEE ALSO
% mmds
% cca (in SOM Toolbox)
% sammon (in SOM Toolbox)

%COPYRIGHT NOTICE
%This function is a part of Icasso software library
%Copyright (C) 2003-2005 Johan Himberg
%
%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 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.

% ver 1.21 030305 johan

% Set default projection method
if nargin<2|isempty(method),
  method='cca';
end

% Check the method
method=lower(method);
switch method
 case {'sammon','cca','mmds'}
  ;
 otherwise
  error('Unknown projection.');
end

% We project onto plane
outputDimension=2;

% Set default parameters for proj, methods
switch method
 case 'sammon'
  default={'alpha',0.7,'epochs',100,'s2d','sqrtsim2dis'};
 case 'cca'
  default={'alpha',0.7,'epochs',75,...
	   'radius',max(icassoGet(sR,'M')/20,10),'s2d','sqrtsim2dis'};
 case 'mmds'
  default={'s2d','sqrtsim2dis'};
end

%% Check optional arguments and add defaults
projectionparameters=processvarargin(varargin,default);
num_of_args=length(projectionparameters);

%% check arguments
for i=1:2:num_of_args;
  switch lower(projectionparameters{i})
   case 's2d'
    sim2dis=projectionparameters{i+1};
   case 'epochs'
    epochs=projectionparameters{i+1};
   case 'alpha'
    alpha=projectionparameters{i+1};
   case 'radius'
    CCAradius=projectionparameters{i+1};
   otherwise
    error(['Indentifier ' projectionparameters{i} ' not recognized.']);
  end
end

% Make similarity-to-dissimilarity transformation

D=feval(sim2dis,sR.cluster.similarity);

disp([char(13) 'Projection, using ' upper(method) char(13)]);

switch method
 case 'mmds'
  P=mmds(D);
  P=P(:,1:outputDimension);
 otherwise
  % Start from MMDS
  initialProjection=mmds(D); initialProjection=initialProjection(:,1:2);
  %
  dummy=rand(size(D,1),outputDimension);
  % rand. init projection: set
  % initialProjection=dummy;
  switch method
   case 'sammon' % Use SOM Toolbox Sammon
    P=sammon(dummy,initialProjection,epochs,'steps',alpha,D);
   case 'cca'    % Use SOM Toolbox CCA
    P=cca(dummy,initialProjection,epochs,D,alpha,CCAradius);
  end
end

sR.projection.method=method;
sR.projection.parameters=projectionparameters;
sR.projection.coordinates=P;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% POP_SELECT - given an input EEG dataset structure, output a new EEG data structure
%                retaining and/or excluding specified time/latency, data point, channel,
%                and/or epoch range(s).
% Usage:
%   >> OUTEEG = pop_select(INEEG, 'key1', value1, 'key2', value2 ...);
%
% Graphic interface:
%   "Time range" - [edit box] RETAIN only the indicated epoch latency or continuous data
%                  time range: [low high] in ms, inclusive. For continuous data, several
%                  time ranges may be specified, separated by semicolons.
%                  Example: "5 10; 12 EEG.xmax" will retain the indicated
%                  stretches of continuous data, and remove data portions outside
%                  the indicated ranges, e.g. from 0 s to 5 s and from 10 s to 12 s.
%                  Command line equivalent: 'time' (or 'notime' - see below)
%   "Time range" - [checkbox] EXCLUDE the indicated latency range(s) from the data.
%                  For epoched data, it is not possible to remove a range of latencies
%                  from the middle of the epoch, so either the low and/or the high values
%                  in the specified latency range (see above) must be at an epoch boundary
%                  (EEG.xmin, EEGxmax).  Command line equivalent: [if checked] 'notime'
%   "Point range" - [edit box] RETAIN the indicated data point range(s).
%                  Same options as for the "Time range" features (above).
%                  Command line equivalent: 'point' (or 'nopoint' - see below).
%   "Point range" - [checkbox] EXCLUDE the indicated point range(s).
%                  Command line equivalent: [if checked] 'nopoint'
%   "Epoch range" - [edit box] RETAIN the indicated data epoch indices in the dataset.
%                  This checkbox is only visible for epoched datasets.
%                  Command line equivalent: 'trial' (or 'notrial' - see below)
%   "Epoch range" - [checkbox] EXCLUDE the specified data epochs.
%                   Command line equivalent: [if checked] 'notrial'
%   "Channel range" - [edit box] RETAIN the indicated vector of data channels
%                  Command line equivalent: 'channel' (or 'nochannel' - see below)
%   "Channel range" - [checkbox] EXCLUDE the indicated channels.
%                  Command line equivalent: [if checked] 'nochannel'
%   "..." - [button] select channels by name.
%   "Scroll dataset" - [button] call the EEGPLOT function to scroll the
%                  channel activities in a new window for visual inspection.
%                  Commandline equivalent: EEGPLOT - see its help for details.
% Inputs:
%   INEEG         - input EEG dataset structure
%
% Optional inputs
%   'time'        - [min max] in seconds. Epoch latency or continuous data time range
%                   to retain in the new dataset, (Note: not ms, as in the GUI text entry
%                   above). For continuous data (only), several time ranges can be specified,
%                   separated by semicolons. Example: "5 10; 12 EEG.xmax" will retain
%                   the indicated times ranges, removing data  outside the indicated ranges
%                   e.g. here from 0 to 5 s and from 10 s to 12 s. (See also, 'notime')
%   'rmtime'      - [min max] in seconds. Epoch latency or continuous dataset time range
%                   to exclude from the new dataset. For continuous data, may be
%                   [min1 max1; min2 max2; ...] to exclude several time ranges. For epoched
%                   data, the latency range must include an epoch boundary, as latency
%                   ranges in the middle of epochs cannot be removed from epoched data.
%   'point'       - [min max] epoch or continuous data point range to retain in the new
%                   dataset. For continuous datasets, this may be [min1 max1; min2 max2; ...]
%                   to retain several point ranges. (Notes: If both 'point'/'nopoint' and
%                   'time' | 'notime' are specified, the 'point' limit values take precedence.
%                   The 'point' argument was originally a point vector, now deprecated).
%   'rmpoint'     - [min max] epoch or continuous data point range to exclude in the new dataset.
%                   For epoched data, the point range must include either the first (0)
%                   or the last point (EEG.pnts), as a central point range cannot be removed.
%   'trial'       - [integer array] array of trial indices to retain in the new dataset
%   'rmtrial'     - [integer array] array of trial indices to exclude from the new dataset
%   'sorttrial'   - ['on'|'off'] sort trial indices before extracting them (default: 'on').
%   'checkchans'  - ['on'|'off'] check that channels are present before
%                   rejecting them (default: 'on')
%   'channel'     - vector of channel indices to retain in the new
%                   dataset. Can also be a cell array of channel names.
%   'rmchannel'   - vector of channel indices to exclude from the new
%                   dataset. Can also be a cell array of channel names.
%   'chantype'    - [string|cell] list of channel types to keep
%   'rmchantype'  - [string|cell] list of channel types to remove
%
% Outputs:
%   OUTEEG        - new EEG dataset structure
%
% Note: This function performs a conjunction (AND) of all its optional inputs.
%       Using negative counterparts of all options, any logical combination is
%       possible. Legacy input 'notrial', 'notime', 'nochannel', 'nopoint'
%       are still supported.
%
% Author: Arnaud Delorme, CNL/Salk Institute, 2001; SCCN/INC/UCSD, 2002-
%
% see also: EEGLAB

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. 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 HOLDER 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.

% 01-25-02 reformated help & license -ad
% 01-26-02 changed the format for events and trial conditions -ad
% 02-04-02 changed display format and allow for negation of inputs -ad
% 02-17-02 removed the event removal -ad
% 03-17-02 added channel info subsets selection -ad
% 03-21-02 added event latency recalculation -ad

function [EEG, com] = pop_select( EEG, varargin)

com = '';
if nargin < 1
    help pop_select;
    return;
end
if isempty(EEG(1).data)
    disp('Pop_select error: cannot process empty dataset'); return;
end

if nargin < 2
   geometry = { [1 1 1] [1 1 0.25 0.23 0.51] [1 1 0.25 0.23 0.51] [1 1 0.25 0.23 0.51] [1 1 0.25 0.23 0.51] ...
           [1 1 0.25 0.23 0.51] [1] [1 1 1]};
   enabletype = ~isempty(EEG(1).chanlocs) && isfield(EEG(1).chanlocs, 'type') && ~isempty(EEG(1).chanlocs(1).type);
   uilist = { ...
         { 'Style', 'text', 'string', 'Select data in:', 'fontweight', 'bold'  }, ...
         { 'Style', 'text', 'string', 'Input desired range', 'fontweight', 'bold'  }, ...
         { 'Style', 'text', 'string', 'on->remove these', 'fontweight', 'bold'  }, ...
         { 'Style', 'text', 'string', 'Time range [min max] (s)', 'fontangle', fastif(length(EEG)>1, 'italic', 'normal') }, ...
         { 'Style', 'edit', 'string', '', 'enable', 'on' }, ...
         { }, { 'Style', 'checkbox', 'string', '    ', 'enable', 'on' },{ }, ...
         ...
         { 'Style', 'text', 'string', 'Point range (ex: [1 10])', 'fontangle', fastif(length(EEG)>1, 'italic', 'normal') }, ...
         { 'Style', 'edit', 'string', '', 'enable', 'on' }, ...
         { }, { 'Style', 'checkbox', 'string', '    ', 'enable', 'on' },{ }, ...
         ...
         { 'Style', 'text', 'string', 'Epoch range (ex: 3:2:10)', 'fontangle', fastif(length(EEG)>1, 'italic', 'normal') }, ...
         { 'Style', 'edit', 'string', '', 'enable', 'on' }, ...
         { }, { 'Style', 'checkbox', 'string', '    ', 'enable', 'on' },{ }, ...
         ...
         { 'Style', 'text', 'string', 'Channel(s)' }, ...
         { 'Style', 'edit', 'string', '', 'tag', 'chans' }, ...
         { }, { 'Style', 'checkbox', 'string', '    ' }, ...
         { 'style' 'pushbutton' 'string'  '...', 'enable' fastif(isempty(EEG(1).chanlocs), 'off', 'on') ...
           'callback' 'pop_chansel(get(gcbf, ''userdata''), ''field'', ''labels'',   ''handle'', findobj(''parent'', gcbf, ''tag'', ''chans''));' }, ...
           ...
         { 'Style', 'text', 'string', 'Channel type(s)' }, ...
         { 'Style', 'edit', 'string', '', 'tag', 'chantype' }, ...
         { }, { 'Style', 'checkbox', 'string', '    ' }, ...
         { 'style' 'pushbutton' 'string'  '...', 'enable' fastif(enabletype, 'on', 'off') ...
           'callback'  'pop_chansel(get(gcbf, ''userdata''), ''field'', ''type'',   ''handle'', findobj(''parent'', gcbf, ''tag'', ''chantype''));' }, ...
         ...
           { }, { }, { 'Style', 'pushbutton', 'string', 'Scroll dataset', 'enable', fastif(length(EEG)>1, 'off', 'on'), 'callback', ...
                          'eegplot(EEG.data, ''srate'', EEG.srate, ''winlength'', 5, ''limits'', [EEG.xmin EEG.xmax]*1000, ''position'', [100 300 800 500], ''xgrid'', ''off'', ''eloc_file'', EEG.chanlocs);' } {}};
%           'callback' 'tmplabels = get(gcbf, ''userdata''); [~, tmpvalchan] = pop_chansel(tmplabels, ''withindex'', ''on''); set(findobj(gcbf, ''tag'', ''chans''), ''string'',tmpvalchan); clear tmplabels tmpvalchan' }, ...
   chanlocs = eeg_mergelocs(EEG.chanlocs);
   results = inputgui( 'geometry', geometry, 'uilist', uilist, 'helpcom', 'pophelp(''pop_select'');', 'title', 'Select data -- pop_select()', 'userdata', chanlocs );
   if isempty(results), return; end

   % decode inputs
   % -------------
   args = {};
   if ~isempty( results{1} )
       if ~results{2}, args = { args{:}, 'time', eval( [ '[' results{1} ']' ] ) };
       else            args = { args{:}, 'rmtime', eval( [ '[' results{1} ']' ] ) }; end
   end

   if ~isempty( results{3} )
       if ~results{4}, args = { args{:}, 'point', eval( [ '[' results{3} ']' ] ) };
       else            args = { args{:}, 'rmpoint', eval( [ '[' results{3} ']' ] ) }; end
   end

   if ~isempty( results{5} )
       if ~results{6}, args = { args{:}, 'trial', eval( [ '[' results{5} ']' ] ) };
       else            args = { args{:}, 'rmtrial', eval( [ '[' results{5} ']' ] ) }; end
   end

   if ~isempty( results{7} )
       [ chaninds, chanlist ] = eeg_decodechan(chanlocs, results{7});
       if isempty(chanlist)
           if length(EEG) > 1 && length(unique([EEG.nbchan])) > 1
               error([ 'Cannot use channel indices when processing multiple datasets' 10 ...
                   'with some channels already removed' ])
           end
           chanlist = chaninds;
       end
       if ~results{8}, args = { args{:}, 'channel'  , chanlist };
       else            args = { args{:}, 'rmchannel', chanlist }; end
   end

   if ~isempty( results{9} )
       [ ~, chantypes ] = eeg_decodechan(chanlocs, results{9}, 'type');
       if ~results{10}, args = { args{:}, 'chantype'  , unique(chantypes) };
       else             args = { args{:}, 'rmchantype', unique(chantypes) }; end
   end

else
    args = varargin;
end

% process multiple datasets
% -------------------------
if length(EEG) > 1
    if nargin < 2
        [ EEG, com ] = eeg_eval( 'pop_select', EEG, 'warning', 'on', 'params', args);
    else
        [ EEG, com ] = eeg_eval( 'pop_select', EEG, 'warning', 'off', 'params',args);
    end
    return;
end

%----------------------------AMICA---------------------------------
if isfield(EEG.etc,'amica') && isfield(EEG.etc.amica,'prob_added')
    for index = 1:2:length(args)
       if strcmpi(args{index}, 'channel')
           args{index+1} = [ args{index+1} EEG.nbchan-(0:2*EEG.etc.amica.num_models-1)];
       end
    end
end
%--------------------------------------------------------------------

g = finputcheck(args, { 'time'    'real'      []         []; ...
                        'notime'  'real'      []         []; ...
                        'rmtime'  'real'      []         []; ...
                        'trial'   'integer'   []         [1:EEG.trials]; ...
                        'notrial' 'integer'   []         []; ...
                        'rmtrial' 'integer'   []         []; ...
                        'point'   'integer'   []         []; ...
                        'nopoint' 'integer'   []         []; ...
                        'rmpoint' 'integer'   []         []; ...
                        'channel'   { 'integer','cell' }  []   [];
                        'nochannel' { 'integer','cell' }   []  [];
                        'rmchannel' { 'integer','cell' }   []  [];
                        'chantype'    { 'string','cell' }    []  {};
                        'rmchantype'  { 'string','cell' }    []  {};
                        'trialcond'   'integer'   []         []; ...
                        'notrialcond' 'integer'   []         []; ...
                        'sort'        'integer'   []         []; ...
                        'sorttrial'   'string'    { 'on','off' } 'on' }, 'pop_select');
if ischar(g), error(g); end
if ~isempty(g.sort)
    if g.sort, g.sorttrial = 'on';
    else       g.sorttrial = 'off';
    end
end
if strcmpi(g.sorttrial, 'on')
    g.trial = sort(setdiff( g.trial, g.notrial ));
    if isempty(g.trial), error('Error: dataset %s is empty',EEG.filename); end
else
    g.trial(ismember(g.trial,g.notrial)) = [];
    % still warn about & remove duplicate trials (may be removed in the future)
    [p,q] = unique_bc(g.trial);
    if length(p) ~= length(g.trial)
        disp('Warning: trial selection contained duplicated elements, which were removed.');
    end
    g.trial = g.trial(sort(q));
end

% rename parameters
% -----------------
if ~isempty(g.rmtime)    g.notime    = g.rmtime; end
if ~isempty(g.rmpoint)   g.nopoint   = g.rmpoint; end
if ~isempty(g.rmtrial)   g.notrial   = g.rmtrial; end
if ~isempty(g.rmchannel) g.nochannel = g.rmchannel; end

% decode channels
% ---------------
if ~isempty(g.channel) || ~isempty(g.nochannel)
    % find channels by name
    if ~isempty(g.channel)
        if ~isempty(g.chantype) || ~isempty(g.rmchantype)
            error('You can select channels by name or by type but not both');
        end
        inds = eeg_decodechan(EEG, g.channel, 'labels', true);
        chanFlag = zeros(1, EEG.nbchan);
        chanFlag(inds) = 1;
    else
        chanFlag = ones(1, EEG.nbchan);
    end
    if ~isempty(g.nochannel)
        if ~isempty(g.chantype) || ~isempty(g.rmchantype)
            error('You can select channels by name or by type but not both');
        end
        inds = eeg_decodechan(EEG, g.nochannel, 'labels', true);
        chanFlag(inds) = 0;
    end
else
    % find channels by type
    if ~isempty(g.chantype)
        inds = eeg_decodechan(EEG, g.chantype, 'type', true);
        chanFlag = zeros(1, EEG.nbchan);
        chanFlag(inds) = 1;
    else
        chanFlag = ones(1, EEG.nbchan);
    end
    if ~isempty(g.rmchantype)
        inds = eeg_decodechan(EEG, g.rmchantype, 'type', true);
        chanFlag(inds) = 0;
    end
end
g.channel = find(chanFlag);

% time selection
% --------------
if ~isempty(g.time) && (g.time(1) < EEG.xmin*1000) && (g.time(2) > EEG.xmax*1000)
   error('Wrong time range');
end
if min(g.trial) < 1 || max( g.trial ) > EEG.trials
   error('Wrong trial range');
end
if size(g.point,2) > 2
    g.point = [g.point(1) g.point(end)];
    disp('Warning: vector format for point range is deprecated');
end
if size(g.nopoint,2) > 2
    g.nopoint = [g.nopoint(1) g.nopoint(end)];
    disp('Warning: vector format for point range is deprecated');
end
if ~isempty( g.point )
    g.time = zeros(size(g.point));
    for index = 1:length(g.point(:))
        g.time(index) = eeg_point2lat(g.point(index), 1, EEG.srate, [EEG.xmin EEG.xmax]);
    end
    g.notime = [];
end
if ~isempty( g.nopoint )
    g.notime = zeros(size(g.nopoint));
    for index = 1:length(g.nopoint(:))
        g.notime(index) = eeg_point2lat(g.nopoint(index), 1, EEG.srate, [EEG.xmin EEG.xmax]);
    end
    g.time = [];
end
if ~isempty( g.notime )
    if size(g.notime,2) ~= 2
        error('Time/point range must contain 2 columns exactly');
    end
    if g.notime(2) == EEG.xmax
        g.time = [EEG.xmin g.notime(1)];
    else
        if g.notime(1) == EEG.xmin
            g.time = [g.notime(2) EEG.xmax];
        elseif EEG.trials > 1
            error('Wrong notime range. Remember that it is not possible to remove a slice of time for data epochs.');
        end
    end
    if g.notime(end) > EEG.xmax, g.notime(end) = EEG.xmax; end
    if g.notime(1)   < EEG.xmin, g.notime(1)   = EEG.xmin; end
    if floor(max(g.notime(:))) > EEG.xmax
        error('Time/point range exceed upper data limits');
    end
    if min(g.notime(:)) < EEG.xmin
        error('Time/point range exceed lower data limits');
    end
end
if ~isempty(g.time)
    if size(g.time,2) ~= 2
        error('Time/point range must contain 2 columns exactly');
    end
    for index = 1:length(g.time)
        if g.time(index) > EEG.xmax
            g.time(index) = EEG.xmax;
            disp('Upper time limits exceed data, corrected');
        elseif g.time(index) < EEG.xmin
            g.time(index) = EEG.xmin;
            disp('Lower time limits exceed data, corrected');
        end
    end
end

% select trial values
%--------------------
if ~isempty(g.trialcond)
   try
        tt = struct( g.trialcond{:} ); catch
        error('Trial conditions format error');
   end
   ttfields = fieldnames (tt);
   for index = 1:length(ttfields)
        if ~isfield( EEG.epoch, ttfields{index} )
            error([ ttfields{index} 'is not a field of EEG.epoch' ]);
        end
        tmpepoch = EEG.epoch;
	    eval( [ 'Itriallow  = find( [ tmpepoch(:).' ttfields{index} ' ] >= tt.' ttfields{index} '(1) );' ] );
	    eval( [ 'Itrialhigh = find( [ tmpepoch(:).' ttfields{index} ' ] <= tt.' ttfields{index} '(end) );' ] );
	    Itrialtmp = intersect_bc(Itriallow, Itrialhigh);
	    g.trial = intersect_bc( g.trial(:)', Itrialtmp(:)');
   end
end

if isempty(g.trial)
   error('Empty dataset, no trial');
end
if length(g.trial) ~= EEG.trials
	fprintf('Removing %d trial(s)...\n', EEG.trials - length(g.trial));
end
if length(g.channel) ~= EEG.nbchan
	fprintf('Removing %d channel(s)...\n', EEG.nbchan - length(g.channel));
end

try
    % For AMICA probabilities...
    %-----------------------------------------------------
    if isfield(EEG.etc, 'amica') && ~isempty(EEG.etc.amica) && isfield(EEG.etc.amica, 'v_smooth') && ~isempty(EEG.etc.amica.v_smooth) && ~isfield(EEG.etc.amica,'prob_added')
        if isfield(EEG.etc.amica, 'num_models') && ~isempty(EEG.etc.amica.num_models)
            if size(EEG.data,2) == size(EEG.etc.amica.v_smooth,2) && size(EEG.data,3) == size(EEG.etc.amica.v_smooth,3) && size(EEG.etc.amica.v_smooth,1) == EEG.etc.amica.num_models

                EEG = eeg_formatamica(EEG);

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

                [EEG, com] = pop_select(EEG,args{:});

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

                EEG = eeg_reformatamica(EEG);
                EEG = eeg_checkamica(EEG);
                return;
            else
                disp('AMICA probabilities not compatible with size of data, probabilities cannot be rejected')

                disp('Resuming rejection...')
            end
        end

    end
    % ------------------------------------------------------
catch
    warnmsg = strcat('your dataset contains amica information, but the amica plugin is not installed.  Continuing and ignoring amica information.');
    warning(warnmsg)
end

% recompute latency and epoch number for events
% ---------------------------------------------
if length(g.trial) ~= EEG.trials && ~isempty(EEG.event)
    if ~isfield(EEG.event, 'epoch')
        disp('Pop_epoch warning: bad event format with epoch dataset, removing events');
        EEG.event = [];
    else
        if isfield(EEG.event, 'epoch')
            keepevent = [];
            for indexevent = 1:length(EEG.event)
                newindex = find( EEG.event(indexevent).epoch == g.trial );% For AMICA probabilities...
                %-----------------------------------------------------
                try
                    if isfield(EEG.etc, 'amica') && ~isempty(EEG.etc.amica) && isfield(EEG.etc.amica, 'v_smooth') && ~isempty(EEG.etc.amica.v_smooth) && ~isfield(EEG.etc.amica,'prob_added')
                        if isfield(EEG.etc.amica, 'num_models') && ~isempty(EEG.etc.amica.num_models)
                            if size(EEG.data,2) == size(EEG.etc.amica.v_smooth,2) && size(EEG.data,3) == size(EEG.etc.amica.v_smooth,3) && size(EEG.etc.amica.v_smooth,1) == EEG.etc.amica.num_models

                                EEG = eeg_formatamica(EEG);

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

                                [EEG, com] = pop_select(EEG,args{:});

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

                                EEG = eeg_reformatamica(EEG);
                                EEG = eeg_checkamica(EEG);
                                return;
                            else
                                disp('AMICA probabilities not compatible with size of data, probabilities cannot be rejected')

                                disp('Resuming rejection...')
                            end
                        end

                    end
                catch
                    warnmsg = strcat('your dataset contains amica information, but the amica plugin is not installed.  Continuing and ignoring amica information.');
                    warning(warnmsg)
                end;
                % ------------------------------------------------------

                if ~isempty(newindex)
                    keepevent = [keepevent indexevent];
                    if isfield(EEG.event, 'latency')
                        EEG.event(indexevent).latency = EEG.event(indexevent).latency - (EEG.event(indexevent).epoch-1)*EEG.pnts + (newindex-1)*EEG.pnts;
                    end
                    EEG.event(indexevent).epoch = newindex;
                end
            end
            diffevent = setdiff_bc([1:length(EEG.event)], keepevent);
            if ~isempty(diffevent)
                disp(['Pop_select: removing ' int2str(length(diffevent)) ' unreferenced events']);
                EEG.event(diffevent) = [];
            end
        end
    end
end

% performing removal
% ------------------
if ~isempty(g.time) || ~isempty(g.notime)
    if EEG.trials > 1
        % select new time window
        % ----------------------
        try,   tmpevent = EEG.event;
               tmpeventlatency = [ tmpevent.latency ];
        catch, tmpeventlatency = [];
        end
        alllatencies = 1-(EEG.xmin*EEG.srate); % time 0 point
        alllatencies = linspace( alllatencies, EEG.pnts*(EEG.trials-1)+alllatencies, EEG.trials);
        [EEG.data, tmptime, indices, epochevent]= epoch(EEG.data, alllatencies, ...
                                                     [g.time(1) g.time(2)]*EEG.srate, 'allevents', tmpeventlatency);
        tmptime = tmptime/EEG.srate;
        if g.time(1) ~= tmptime(1) && g.time(2)-1/EEG.srate ~= tmptime(2)
            fprintf('pop_select(): time limits have been adjusted to [%3.3f %3.3f] to fit data points limits\n', tmptime(1), tmptime(2)+1/EEG.srate);
        end
        EEG.xmin = tmptime(1);
        EEG.xmax = tmptime(2);
        EEG.pnts = size(EEG.data,2);
        alllatencies = alllatencies(indices);

        % modify the event structure accordingly (latencies and add epoch field)
        % ----------------------------------------------------------------------
        allevents = [];
        newevent = [];
        count = 1;
        if ~isempty(epochevent)
            newevent = EEG.event(1);
            for index=1:EEG.trials
                for indexevent = epochevent{index}
                    newevent(count)         = EEG.event(indexevent);
                    newevent(count).epoch   = index;
                    newevent(count).latency = newevent(count).latency - alllatencies(index) - tmptime(1)*EEG.srate + 1 + EEG.pnts*(index-1);
                    count = count + 1;
                end
            end
        end
        EEG.event = newevent;

        % erase event-related fields from the epochs
        % ------------------------------------------
        if ~isempty(EEG.epoch)
            fn = fieldnames(EEG.epoch);
            EEG.epoch = rmfield(EEG.epoch,{fn{strmatch('event',fn)}});
        end
    else
        if isempty(g.notime)
            if length(g.time) == 2 && EEG.xmin < 0
                disp('Warning: negative minimum time; unchanged to ensure correct latency of initial boundary event');
            end
            g.notime = g.time';
            g.notime = g.notime(:);
            if g.notime(1) ~= 0, g.notime = [EEG.xmin g.notime(:)'];
            else                 g.notime = [g.notime(2:end)'];
            end
            if g.time(end) == EEG.xmax, g.notime(end) = [];
            else                        g.notime(end+1) = EEG.xmax;
            end

            for index = 1:length(g.notime)
                if g.notime(index) ~= 0  && g.notime(index) ~= EEG.xmax
                    if mod(index,2), g.notime(index) = g.notime(index) + 1/EEG.srate;
                    else             g.notime(index) = g.notime(index) - 1/EEG.srate;
                    end
                end
            end
            g.notime = reshape(g.notime, 2, length(g.notime)/2)';
        end

        nbtimes = length(g.notime(:));
        [points,flag] = eeg_lat2point(g.notime(:)', ones(1,nbtimes), EEG.srate, [EEG.xmin EEG.xmax]);
        points = reshape(points, size(g.notime));

        % fixing if last region is the same
        if flag
            if ~isempty(find((points(end,1)-points(end,2))== 0)), points(end,:) = []; end
        end

        EEG = eeg_eegrej(EEG, points);
    end
end

% performing removal
% ------------------
if ~isequal(g.channel,1:size(EEG.data,1)) || ~isequal(g.trial,1:size(EEG.data,3))
    eeglab_options;
    if option_memmapdata
        % this code below is preferred for memory mapped files
        diff1 = setdiff_bc([1:size(EEG.data,1)], g.channel);
        diff2 = setdiff_bc([1:size(EEG.data,3)], g.trial);
        if ~isempty(diff1)
            EEG.data(diff1, :, :) = [];
        end
        if ~isempty(diff2)
            EEG.data(:, :, diff2) = [];
        end
    else
        EEG.data  = EEG.data(g.channel, :, g.trial);
    end
end
if ~isempty(EEG.icaact), EEG.icaact = EEG.icaact(:,:,g.trial); end
if ~isempty(EEG.chanlocs)
    if ~isfield(EEG.chaninfo, 'removedchans')
        EEG.chaninfo.removedchans = [];
    end
    try
        diff1 = setdiff_bc([1:EEG.nbchan], g.channel);
        fields = fieldnames(EEG.chanlocs);
        for iChan = diff1(:)'
            EEG.chaninfo.removedchans(end+1).(fields{1}) = EEG.chanlocs(iChan).(fields{1});
            for iField = 1:length(fields)
                EEG.chaninfo.removedchans(end).(fields{iField}) = EEG.chanlocs(iChan).(fields{iField});
            end
        end
    catch
        disp('There was an issue storing removed channels in pop_select');
    end
    EEG.chanlocs = EEG.chanlocs(g.channel);
end
EEG.trials    = length(g.trial);
EEG.pnts      = size(EEG.data,2);
EEG.nbchan    = length(g.channel);
if ~isempty(EEG.epoch)
   EEG.epoch = EEG.epoch( g.trial );
end
if ~isempty(EEG.specdata)
	if length(g.point) == EEG.pnts
   		EEG.specdata = EEG.specdata(g.channel, :, g.trial);
   	else
   		EEG.specdata = [];
   		fprintf('Warning: spectral data were removed because of the change in the number of points\n');
    end
end

% ica specific
% ------------
if ~isempty(EEG.icachansind)

    rmchans = setdiff_bc( EEG.icachansind, g.channel ); % channels to remove

    % channel sub-indices
    % -------------------
    icachans = 1:length(EEG.icachansind);
    for index = length(rmchans):-1:1
        chanind           = find(EEG.icachansind == rmchans(index));
        icachans(chanind) = [];
    end

    % new channels indices
    % --------------------
    count   = 1;
    newinds = [];
    for index = 1:length(g.channel)
        if any(EEG.icachansind == g.channel(index))
            newinds(count) = index;
            count          = count+1;
        end
    end
    EEG.icachansind = newinds;

else
    icachans = 1:size(EEG.icasphere,2);
end

if ~isempty(EEG.icawinv)
    flag_rmchan = (length(icachans) ~= size(EEG.icawinv,1));
    if  isempty(EEG.icaweights) || flag_rmchan
        EEG.icawinv    = EEG.icawinv(icachans,:);
        EEG.icaweights = pinv(EEG.icawinv);
        EEG.icasphere  = eye(size(EEG.icaweights,2));
    end
end
if ~isempty(EEG.specicaact)
    if length(g.point) == EEG.pnts
        EEG.specicaact = EEG.specicaact(icachans, :, g.trial);
    else
        EEG.specicaact = [];
        fprintf('Warning: spectral ICA data were removed because of the change in the number of points\n');
    end
end

% check if only one epoch
% -----------------------
if EEG.trials == 1
    if isfield(EEG.event, 'epoch')
        EEG.event = rmfield(EEG.event, 'epoch');
    end
    EEG.epoch = [];
end
if isfield(EEG.reject, 'gcompreject') && isequal(g.channel,1:size(EEG.data,1))
    tmpgcompreject = EEG.reject.gcompreject;
    EEG.reject = [];
    EEG.reject.gcompreject = tmpgcompreject;
else
    EEG.reject = [];
end
EEG.stats  = [];
EEG.reject.rejmanual = [];
% for stats, can adapt remove the selected trials and electrodes
% in the future to gain time -----------------------------------
EEG.stats.jp = [];
EEG = eeg_checkset(EEG, 'eventconsistency');

% generate command
% ----------------
if nargout > 1
    com = sprintf('EEG = pop_select( EEG, %s);', vararg2str(args));
end

return;

% ********* OLD, do not remove any event any more
% ********* in the future maybe do a pack event to remove events not in the time range of any epoch

if ~isempty(EEG.event)
    % go to array format if necessary
    if isstruct(EEG.event), format = 'struct';
    else                     format = 'array';
    end
    switch format, case 'struct', EEG = eventsformat(EEG, 'array'); end

    % keep only events related to the selected trials
    Indexes = [];
    Ievent  = [];
    for index = 1:length( g.trial )
        currentevents = find( EEG.event(:,2) == g.trial(index));
        Indexes = [ Indexes ones(1, length(currentevents))*index ];
        Ievent  = union_bc( Ievent, currentevents );
    end
    EEG.event = EEG.event( Ievent,: );
    EEG.event(:,2) = Indexes(:);

    switch format, case 'struct', EEG = eventsformat(EEG, 'struct'); end
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function X = lmvuEmbed(Y,dims,k,nr_landmark)

% LMVUEMBED Embed data set with landmark MVU
% FORMAT
% DESC Embed data set with landmark version of Weinberg et al.'s
% maximum variance unfolding algorithm.
% ARG Y : Data
% ARG dims : Dimensionality of Embedding (default = 2)
% ARG k : Number of Neighbours in Proximity Graph (default = 7)
% ARG nr_landmark : Number of landmark Points
% RETURN X : embedding
%
% SEEALSO : ppcaEmbed, lleEmbed, mvuEmbed
%
% COPYRIGHT : Carl Henrik Ek, Neil D. Lawrence, 2007

% MLTOOLS

if(nargin<4)
  nr_landmark = 30;
  if(nargin<3)
    k = 7;
    if(nargin<2)
      dims = 2;
      if(nargin<1)
	error('To Few Arguments');
      end
    end
  end
end

if(any(any(isnan(Y))))
  error('Cannot Initialise GPLVM using lmvu when missing data is present.');
end

X = lmvu(distance(Y'),nr_landmark,k);

X = X(1:1:dims,:)';

return

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
if ~exist('save_name'),
    save_name = 'Calib_Results';
end;

fprintf(1,'Generating the matlab script file %s.m containing the intrinsic and extrinsic parameters...\n',save_name)

fid = fopen([ save_name '.m'],'wt');

fprintf(fid,'%% Intrinsic and Extrinsic Camera Parameters\n');
fprintf(fid,'%%\n');
fprintf(fid,'%% This script file can be directly excecuted under Matlab to recover the camera intrinsic and extrinsic parameters.\n');
fprintf(fid,'%% IMPORTANT: This file contains neither the structure of the calibration objects nor the image coordinates of the calibration points.\n');
fprintf(fid,'%%            All those complementary variables are saved in the complete matlab data file Calib_Results.mat.\n');
fprintf(fid,'%% For more information regarding the calibration model visit http://www.vision.caltech.edu/bouguetj/calib_doc/\n');
fprintf(fid,'\n\n');
fprintf(fid,'%%-- Focal length:\n');
fprintf(fid,'fc = [ %5.15f ; %5.15f ];\n',fc);
fprintf(fid,'\n');
fprintf(fid,'%%-- Principal point:\n');
fprintf(fid,'cc = [ %5.15f ; %5.15f ];\n',cc);
fprintf(fid,'\n');
fprintf(fid,'%%-- Skew coefficient:\n');
fprintf(fid,'alpha_c = %5.15f;\n',alpha_c);
fprintf(fid,'\n');
fprintf(fid,'%%-- Distortion coefficients:\n');
fprintf(fid,'kc = [ %5.15f ; %5.15f ; %5.15f ; %5.15f ; %5.15f ];\n',kc);
fprintf(fid,'\n');
fprintf(fid,'%%-- Focal length uncertainty:\n');
fprintf(fid,'fc_error = [ %5.15f ; %5.15f ];\n',fc_error);
fprintf(fid,'\n');
fprintf(fid,'%%-- Principal point uncertainty:\n');
fprintf(fid,'cc_error = [ %5.15f ; %5.15f ];\n',cc_error);
fprintf(fid,'\n');
fprintf(fid,'%%-- Skew coefficient uncertainty:\n');
fprintf(fid,'alpha_c_error = %5.15f;\n',alpha_c_error);
fprintf(fid,'\n');
fprintf(fid,'%%-- Distortion coefficients uncertainty:\n');
fprintf(fid,'kc_error = [ %5.15f ; %5.15f ; %5.15f ; %5.15f ; %5.15f ];\n',kc_error);
fprintf(fid,'\n');
fprintf(fid,'%%-- Image size:\n');
fprintf(fid,'nx = %d;\n',nx);
fprintf(fid,'ny = %d;\n',ny);
fprintf(fid,'\n');
fprintf(fid,'\n');
fprintf(fid,'%%-- Various other variables (may be ignored if you do not use the Matlab Calibration Toolbox):\n');
fprintf(fid,'%%-- Those variables are used to control which intrinsic parameters should be optimized\n');
fprintf(fid,'\n');
fprintf(fid,'n_ima = %d;\t\t\t\t\t\t%% Number of calibration images\n',n_ima);
fprintf(fid,'est_fc = [ %d ; %d ];\t\t\t\t\t%% Estimation indicator of the two focal variables\n',est_fc);
fprintf(fid,'est_aspect_ratio = %d;\t\t\t\t%% Estimation indicator of the aspect ratio fc(2)/fc(1)\n',est_aspect_ratio);
fprintf(fid,'center_optim = %d;\t\t\t\t\t%% Estimation indicator of the principal point\n',center_optim);
fprintf(fid,'est_alpha = %d;\t\t\t\t\t\t%% Estimation indicator of the skew coefficient\n',est_alpha);
fprintf(fid,'est_dist = [ %d ; %d ; %d ; %d ; %d ];\t%% Estimation indicator of the distortion coefficients\n',est_dist);
fprintf(fid,'\n\n');
fprintf(fid,'%%-- Extrinsic parameters:\n');
fprintf(fid,'%%-- The rotation (omc_kk) and the translation (Tc_kk) vectors for every calibration image and their uncertainties\n');
fprintf(fid,'\n');
for kk = 1:n_ima,
    fprintf(fid,'%%-- Image #%d:\n',kk);
    eval(['omckk = omc_' num2str(kk) ';']);
    eval(['Tckk = Tc_' num2str(kk) ';']);
    fprintf(fid,'omc_%d = [ %d ; %d ; %d ];\n',kk,omckk);
    fprintf(fid,'Tc_%d  = [ %d ; %d ; %d ];\n',kk,Tckk);
    if (exist(['Tc_error_' num2str(kk)])==1) & (exist(['omc_error_' num2str(kk)])==1),
        eval(['omckk_error = omc_error_' num2str(kk) ';']);
        eval(['Tckk_error = Tc_error_' num2str(kk) ';']);
        fprintf(fid,'omc_error_%d = [ %d ; %d ; %d ];\n',kk,omckk_error);
        fprintf(fid,'Tc_error_%d  = [ %d ; %d ; %d ];\n',kk,Tckk_error);
    end;
    fprintf(fid,'\n');
end;

fclose(fid);

repo: ekfmonoslam


text:
function [elem,bdFlag] = sortelem3(elem,bdFlag)
%% SORTELEM3 sort elem in ascend ordering
%
% [elem,bdFlag] = sortelem3(elem,bdFlag) sorts the elem such that
% elem(t,1)< elem(t,2)< elem(t,3)<elem(t,4). A simple sort(elem,2) cannot
% sort bdFlag.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

%% Step 1: elem(:,4) is the biggest one
[tempvar,idx] = max(elem,[],2);  %#ok<*ASGLU>
elem(idx==1,1:4) = elem(idx==1,[2 4 3 1]);
elem(idx==2,1:4) = elem(idx==2,[3 4 1 2]);
elem(idx==3,1:4) = elem(idx==3,[4 2 1 3]);
if exist('bdFlag','var')
    bdFlag(idx==1,1:4) = bdFlag(idx==1,[2 4 3 1]);
    bdFlag(idx==2,1:4) = bdFlag(idx==2,[3 4 1 2]);
    bdFlag(idx==3,1:4) = bdFlag(idx==3,[4 2 1 3]);
end
%% Step 2: elem(:,1) is the smallest one
[tempvar,idx] = min(elem(:,1:3),[],2);
% elem(idx==1,1:3) = elem(idx==1,[1 2 3]);
elem(idx==2,1:3) = elem(idx==2,[2 3 1]);
elem(idx==3,1:3) = elem(idx==3,[3 1 2]);
if exist('bdFlag','var')
    bdFlag(idx==2,1:3) = bdFlag(idx==2,[2 3 1]);
    bdFlag(idx==3,1:3) = bdFlag(idx==3,[3 1 2]);
end

%% Step 3: sort elem(:,2)<elem(:3)
idx = (elem(:,3) < elem(:,2));
elem(idx,[2 3]) = elem(idx,[3 2]);
if exist('bdFlag','var')
    bdFlag(idx,[2 3]) = bdFlag(idx,[3 2]);
end

%% Output
if ~exist('bdFlag','var')
    bdFlag = [];
end

repo: Aero-matlab


text:
function [varargout] = wavefilter(wname, type)
%WAVEFILTER Create wavelet decomposition and reconstruction filters.
%   [VARARGOUT] = WAVEFILTER(WNAME, TYPE) returns the decomposition
%   and/or reconstruction filters used in the computation of the
%   forward and inverse FWT (fast wavelet transform).
%
%   EXAMPLES:
%     [ld, hd, lr, hr] = wavefilter('haar') Get the low and highpass
%                                           decomposition (ld, hd)
%                                           and reconstruction
%                                           (lr, hr) filters for
%                                           wavelet 'haar'.
%     [ld, hd] = wavefilter('haar','d')     Get decomposition filters
%                                           ld and hd.
%     [lr, hr] = wavefilter('haar','r')     Get reconstruction
%                                           filters lr and hr.
%
%   INPUTS:
%     WNAME             Wavelet Name
%     ---------------------------------------------------------
%     'haar' or 'db1'   Haar
%     'db4'             4th order Daubechies
%     'sym4'            4th order Symlets
%     'bior6.8'         Cohen-Daubechies-Feauveau biorthogonal
%     'jpeg9.7'         Antonini-Barlaud-Mathieu-Daubechies
%
%     TYPE              Filter Type
%     ---------------------------------------------------------
%     'd'               Decomposition filters
%     'r'               Reconstruction filters
%
%   See also WAVEFAST and WAVEBACK.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.5 $  $Date: 2003/10/13 01:09:39 $

% Check the input and output arguments.
error(nargchk(1, 2, nargin));

if (nargin == 1 & nargout ~= 4) | (nargin == 2 & nargout ~= 2)
   error('Invalid number of output arguments.');
end

if nargin == 1 & ~ischar(wname)
   error('WNAME must be a string.');
end

if nargin == 2 & ~ischar(type)
   error('TYPE must be a string.');
end

% Create filters for the requested wavelet.
switch lower(wname)
case {'haar', 'db1'}
   ld = [1 1]/sqrt(2);     hd = [-1 1]/sqrt(2);
   lr = ld;                hr = -hd;

case 'db4'
   ld = [-1.059740178499728e-002 3.288301166698295e-002 ...
         3.084138183598697e-002 -1.870348117188811e-001 ...
         -2.798376941698385e-002 6.308807679295904e-001 ...
         7.148465705525415e-001 2.303778133088552e-001];
   t = (0:7);
   hd = ld;    hd(end:-1:1) = cos(pi * t) .* ld;
   lr = ld;    lr(end:-1:1) = ld;
   hr = cos(pi * t) .* ld;

case 'sym4'
   ld = [-7.576571478927333e-002 -2.963552764599851e-002 ...
         4.976186676320155e-001 8.037387518059161e-001 ...
         2.978577956052774e-001 -9.921954357684722e-002 ...
         -1.260396726203783e-002 3.222310060404270e-002];
   t = (0:7);
   hd = ld;    hd(end:-1:1) = cos(pi * t) .* ld;
   lr = ld;    lr(end:-1:1) = ld;
   hr = cos(pi * t) .* ld;

case 'bior6.8'
   ld = [0 1.908831736481291e-003 -1.914286129088767e-003 ...
         -1.699063986760234e-002 1.193456527972926e-002 ...
         4.973290349094079e-002 -7.726317316720414e-002 ...
         -9.405920349573646e-002 4.207962846098268e-001 ...
         8.259229974584023e-001 4.207962846098268e-001 ...
         -9.405920349573646e-002 -7.726317316720414e-002 ...
         4.973290349094079e-002 1.193456527972926e-002 ...
         -1.699063986760234e-002 -1.914286129088767e-003 ...
         1.908831736481291e-003];
   hd = [0 0 0 1.442628250562444e-002 -1.446750489679015e-002 ...
         -7.872200106262882e-002 4.036797903033992e-002 ...
         4.178491091502746e-001 -7.589077294536542e-001 ...
         4.178491091502746e-001 4.036797903033992e-002 ...
         -7.872200106262882e-002 -1.446750489679015e-002 ...
         1.442628250562444e-002 0 0 0 0];
   t = (0:17);
   lr = cos(pi * (t + 1)) .* hd;
   hr = cos(pi * t) .* ld;

case 'jpeg9.7'
   ld = [0 0.02674875741080976 -0.01686411844287495 ...
         -0.07822326652898785 0.2668641184428723 ...
         0.6029490182363579 0.2668641184428723 ...
         -0.07822326652898785 -0.01686411844287495 ...
         0.02674875741080976];
   hd = [0 -0.09127176311424948 0.05754352622849957 ...
         0.5912717631142470 -1.115087052456994 ...
         0.5912717631142470 0.05754352622849957 ...
         -0.09127176311424948 0 0];
   t = (0:9);
   lr = cos(pi * (t + 1)) .* hd;
   hr = cos(pi * t) .* ld;

otherwise
   error('Unrecognizable wavelet name (WNAME).');
end

% Output the requested filters.
if (nargin == 1)
   varargout(1:4) = {ld, hd, lr, hr};
else
   switch lower(type(1))
   case 'd'
      varargout = {ld, hd};
   case 'r'
      varargout = {lr, hr};
   otherwise
      error('Unrecognizable filter TYPE.');
   end
end

repo: weiyanmin


text:
function [p,f] = mrv_spm_powell(p,xi,tolsc,func,varargin)
%
% NOTE: this is a slightly modified version of spm_powell. It is
% apparently not compatible with the original.
%
% Powell optimisation method
% FORMAT [p,f] = mrv_spm_powell(p,xi,tolsc,func,varargin)
% 	p        - Starting parameter values
%	xi       - columns containing directions in which to begin
% 	           searching.
% 	tolsc    - stopping criteria
% 	           - optimisation stops when
% 	             sqrt(sum(((p-p_prev)./tolsc).^2))<1
% 	func     - name of evaluated function
% 	varargin - remaining arguments to func (after p)
%
% 	p        - final parameter estimates
% 	f        - function value at minimum
%
%_______________________________________________________________________
% Method is based on Powell's optimisation method from Numerical Recipes
% in C (Press, Flannery, Teukolsky & Vetterling).
%_______________________________________________________________________
% @(#)spm_powell.m	2.3 John Ashburner 01/09/28

p     = p(:);
f     = feval(func,p,varargin{:});
for iter=1:128,
	fprintf('iteration %d...\n', iter);
	pp   = p;
	fp   = f;
	del  = 0;
	for i=1:length(p),
		ft = f;
		[p,junk,f] = linmin(p,xi(:,i),func,f,tolsc,varargin{:});
		if abs(ft-f) > del,
			del  = abs(ft-f);
			ibig = i;
		end;
	end;
	if sqrt(sum(((p(:)-pp(:))./tolsc(:)).^2))<1, return; end;
	ft = feval(func,2.0*p-pp,varargin{:});
	if ft < f,
		[p,xi(:,ibig),f] = linmin(p,p-pp,func,f,tolsc,varargin{:});
	end;
end;
warning('Too many optimisation iterations');
return;
%_______________________________________________________________________

%_______________________________________________________________________
function [p,pi,f] = linmin(p,pi,func,f,tolsc,varargin)
% Line search for minimum.

global lnm % used in linmineval
lnm      = struct('p',p,'pi',pi,'func',func,'args',[]);
lnm.args = varargin;

linmin_plot('Init', 'Line Minimisation','Function','Parameter Value');
linmin_plot('Set', 0, f);

tol      = 1/sqrt(sum((pi(:)./tolsc(:)).^2));
t        = bracket(f);
[f,pmin] = brents(t,tol);
pi       = pi*pmin;
p        = p + pi;

for i=1:length(p), fprintf('%-8.4g ', p(i)); end;
fprintf('| %.5g\n', f);
linmin_plot('Clear');

return;
%_______________________________________________________________________

%_______________________________________________________________________
function f = linmineval(p)
% Reconstruct parameters and evaluate.

global lnm % defined in linmin
pt = lnm.p+p.*lnm.pi;
f  = feval(lnm.func,pt,lnm.args{:});
linmin_plot('Set',p,f);
return;
%_______________________________________________________________________

%_______________________________________________________________________
function t = bracket(f)
% Bracket the minimum (t(2)) between t(1) and t(2)

gold   = (1+sqrt(5))/2; % Golden ratio

t(1)   = struct('p',0,'f',f);
t(2).p = 1;
t(2).f = linmineval(t(2).p);

% if not better then swap
if t(2).f > t(1).f,
	tmp  = t(1);
	t(1) = t(2);
	t(2) = tmp;
end;

t(3).p = t(2).p + gold*(t(2).p-t(1).p);
t(3).f = linmineval(t(3).p);

while t(2).f > t(3).f,

	% fit a polynomial to t
	tmp = cat(1,t.p)-t(2).p;
	pol = pinv([ones(3,1) tmp tmp.^2])*cat(1,t.f);

	% minimum is when gradient of polynomial is zero
	% sign of pol(3) (the 2nd deriv) should be +ve
	d   = -pol(2)/(2*pol(3)+eps);
	u.p = t(2).p+d;

	ulim = t(2).p+100*(t(3).p-t(2).p);
	if (t(2).p-u.p)*(u.p-t(3).p) > 0.0,
		% u is between t(2) and t(3)
		u.f = linmineval(u.p);
		if u.f < t(3).f,
			% minimum between t(2) and t(3) - done
			t(1) = t(2);
			t(2) = u;
			return;
		elseif u.f > t(2).f,
			% minimum between t(1) and u - done
			t(3) = u;
			return;
		end;
		% try golden search instead
		u.p = t(3).p+gold*(t(3).p-t(2).p);
		u.f = linmineval(u.p);

	elseif (t(3).p-u.p)*(u.p-ulim) > 0.0
		% u is between t(3) and ulim
		u.f = linmineval(u.p);
		if u.f < t(3).f,
			% still no minimum as function is still decreasing
			% t(1) = t(2);
			t(2) = t(3);
			t(3) = u;
			u.p  = t(3).p+gold*(t(3).p-t(2).p);
			u.f  = linmineval(u.p);
		end;

	elseif (u.p-ulim)*(ulim-t(3).p) >= 0.0,
		% gone too far - constrain it
		u.p = ulim;
		u.f = linmineval(u.p);

	else,
		% try golden search instead
		u.p = t(3).p+gold*(t(3).p-t(2).p);
		u.f = linmineval(u.p);
	end;

	% Move all 3 points along
	t(1) = t(2);
	t(2) = t(3);
	t(3) = u;
end;
return;
%_______________________________________________________________________

%_______________________________________________________________________
function [f,p] = brents(t, tol)
% Brent's method for line searching - given that minimum is bracketed

% 1 - golden ratio
Cgold = 1 - (sqrt(5)-1)/2;

% Current and previous displacements
d     = Inf;
pd    = Inf;

% t(1) and t(3) bracket the minimum
if t(1).p>t(3).p,
	brk(1) = t(3).p;
	brk(2) = t(1).p;
else,
	brk(1) = t(1).p;
	brk(2) = t(3).p;
end;

% sort t into best first order
tmp  = t(1);
t(1) = t(2);
t(2) = tmp;
if t(2).f>t(3).f,
	tmp  = t(2);
	t(2) = t(3);
	t(3) = tmp;
end;

for iter=1:128,
	% check stopping criterion
	if abs(t(1).p - 0.5*(brk(1)+brk(2)))+0.5*(brk(2)-brk(1)) <= 2*tol,
		p = t(1).p;
		f = t(1).f;
		return;
	end;

	% keep last two displacents
	ppd = pd;
	pd  = d;

	% fit a polynomial to t
	tmp = cat(1,t.p)-t(1).p;
	pol = pinv([ones(3,1) tmp tmp.^2])*cat(1,t.f);

	% minimum is when gradient of polynomial is zero
	d   = -pol(2)/(2*pol(3)+eps);
	u.p = t(1).p+d;

	% check so that displacement is less than the last but two,
	% that the displaced point is between the brackets
	% and (not sure if it is necessary) that the solution is a minimum
	% rather than a maximum
	eps2 = 2*eps*abs(t(1).p)+eps;
	if abs(d) >= abs(ppd)/2 | u.p <= brk(1)+eps2 | u.p >= brk(2)-eps2 | pol(3)<=0,
		% if criteria are not met, then golden search into the larger part
		if t(1).p >= 0.5*(brk(1)+brk(2)),
			d = Cgold*(brk(1)-t(1).p);
		else,
			d = Cgold*(brk(2)-t(1).p);
		end;
		u.p = t(1).p+d;
	end;

	% FUNCTION EVALUATION
	u.f = linmineval(u.p);

	% Insert the new point into the appropriate position and update
	% the brackets if necessary
	if u.f <= t(1).f,
		if u.p >= t(1).p, brk(1)=t(1).p; else, brk(2)=t(1).p; end;
		t(3) = t(2);
		t(2) = t(1);
		t(1) = u;
	else,
		if u.p < t(1).p, brk(1)=u.p; else, brk(2)=u.p; end;
		if u.f <= t(2).f | t(1).p==t(2).p,
			t(3) = t(2);
			t(2) = u;
		elseif u.f <= t(3).f | t(1).p==t(3).p | t(2).p==t(3).p,
			t(3) = u;
		end;
	end;
end;
warning('Too many iterations in Brents');
return;
%_______________________________________________________________________

%_______________________________________________________________________
function linmin_plot(action,arg1,arg2,arg3,arg4)
% Visual output for line minimisation
global linminplot
%-----------------------------------------------------------------------
if (nargin == 0)
	linmin_plot('Init');
else
	% initialize
	%---------------------------------------------------------------
	if (strcmp(lower(action),'init'))
		if (nargin<4)
			arg3 = 'Function';
			if (nargin<3)
				arg2 = 'Value';
				if (nargin<2)
					arg1 = 'Line minimisation';
				end
			end
		end
		fg = spm_figure('FindWin','Interactive');

		if ~isempty(fg)
			linminplot = struct('pointer',get(fg,'Pointer'),'name',get(fg,'Name'),'ax',[]);
			linmin_plot('Clear');
			set(fg,'Pointer','watch');
			% set(fg,'Name',arg1);
			linminplot.ax = axes('Position', [0.15 0.1 0.8 0.75],...
				'Box', 'on','Parent',fg);
			lab = get(linminplot.ax,'Xlabel');
			set(lab,'string',arg3,'FontSize',10);
			lab = get(linminplot.ax,'Ylabel');
			set(lab,'string',arg2,'FontSize',10);
			lab = get(linminplot.ax,'Title');
			set(lab,'string',arg1);
			line('Xdata',[], 'Ydata',[],...
				'LineWidth',2,'Tag','LinMinPlot','Parent',linminplot.ax,...
				'LineStyle','-','Marker','o');
			drawnow;
		end

	% reset
	%---------------------------------------------------------------
	elseif (strcmp(lower(action),'set'))
		F = spm_figure('FindWin','Interactive');
		br = findobj(F,'Tag','LinMinPlot');
		if (~isempty(br))
			[xd,indx] = sort([get(br,'Xdata') arg1]);
			yd = [get(br,'Ydata') arg2];
			yd = yd(indx);
			set(br,'Ydata',yd,'Xdata',xd);
			drawnow;
		end

	% clear
	%---------------------------------------------------------------
	elseif (strcmp(lower(action),'clear'))
		fg = spm_figure('FindWin','Interactive');
		if isstruct(linminplot),
			if ishandle(linminplot.ax), delete(linminplot.ax); end;
			set(fg,'Pointer',linminplot.pointer);
			set(fg,'Name',linminplot.name);
		end;
		spm_figure('Clear',fg);
		drawnow;
	end;
end
%_______________________________________________________________________

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% TD | CP-APR | PARAFAC/CP decomposition solved by Alternating Poisson Regression
% process_video('TD', 'CP-APR', 'dataset/demo.avi', 'output/demo_CP-APR.avi');

r = 10;
A = double(T);
L = double(cp_apr(T,r));
S = A - L;

author: andrewssobral

repo: lrslibrary


text:
function [soln,eqn,info] = StokesBDM1B(node,elem,bdFlag,pde,option)
%% STOKESRT0 Stokes equations: the lowest order BDM element in 2D.
%
%  [soln,eqn,info] = StokesBDM1B(node,elem,bdFlag,pde,option)
%  uses the lowest order of BDM1 element with bubble function
%  (BDM1B) to approximate velocity u and piecewise
%  constant element (P0)  to approximate pressure p, repectively.
%
%  We solve the following equation:
%       - grad div u + curl rot u + grad p  = f   in \Omega
%                                 - div u   = 0   in \Omega
%                                       u   = g   on \Gamma
%
% Based on a version by Ming Wang. Revised by Lin Zhong. Discussed with Jie
% Zhou and Long Chen. Further clean up by Long Chen. Add MG by Long.

if ~exist('option','var'), option = []; end

%% Mesh and data structure
elemunSort = elem;
[elem,bdFlag] = sortelem(elemunSort,bdFlag);
[elem2edge,edge] = dofedge(elem);
[Clambda,area,elemSign] = curlbasis(node,elem);
elem2dof = [elem size(node,1) + elem2edge]; % P2 elements
edge = double(edge);
elem2edge = double(elem2edge);

%% Assembling Matrix
N = size(node,1); NT = size(elem,1); NE = size(edge,1);
Nw = N+NE+NT; Nu = 2*NE+NT; Np = NT; Ndof = Nu + Np;

t = cputime;
% Mvb: Mass matrix with lumping for vertex: P2 element + bubble
Mvb = getmassmatrixP2(elem2dof,area,'NBB');
invMvb = spdiags(1./diag(Mvb),0,Nw,Nw);
Mv = getmassmatrixP2(elem2dof,area,'NB');

% Meb: Mass matrix for edge: BDM1B element
Meb =  getmassmatvec(elem2edge,area,Clambda,'BDM1B');

% invMt: the inverse of Mass matrix for P0 element
invMt = spdiags(1./area,0,NT,NT);

% B: -divergence operator
RT0B = -icdmat(double(elem2edge),elemSign*[1 -1 1]);% -div for RT0
B = [RT0B sparse(NT,NE+NT)];

% C: curl operator
% curl for RT0
RT0C = icdmat(double(edge),[-1 1]);

% curl in HB (hierat = cputime bases)
HC = blkdiag(RT0C,4*speye(NE,NE),speye(NT,NT));

% -------------   The transform matrix between bases  -------------
% NB                                             HB
% 1) (2\lambda_i -1)\lambda_i +1/9w_b            1) \lambda_i
% 2) 4\lambda_i \lambda_j -4/9w_b         <===>  2) 4\lambda_i\lambda_j
% 3) 27\lambda_1\lambda_2\lambda_3               3) 27\lambda_1\lambda_2\lambda_3
%
% NB = HB*T, where T is the tranfer matrix
%   T = | 1      0    0 |
%       |1/2     1    0 |
%       |1/9    4/9   1 |
N2NE = sparse([edge(:,1);edge(:,2)],repmat(1:NE,1,2),1,N,NE);
N2NT = sparse([elem(:,1);elem(:,2);elem(:,3)], repmat(1:NT,1,3),1,N,NT);
NE2NT = sparse([elem2edge(:,1);elem2edge(:,2);elem2edge(:,3)],repmat(1:NT,1,3),1,NE,NT);
T = [speye(N,N) sparse(N,NE) sparse(N,NT);
    -1/2*N2NE'  speye(NE,NE)   sparse(NE,NT);
    1/9*N2NT' -4/9*NE2NT' speye(NT,NT)];
C = HC*T;

% R: weak rot operator
R = invMvb*C'*Meb;

% Vector Laplacian
A = B'*invMt*B + R'*Mvb*R;

%% Assemble right hand side
localEdge = [2,3; 1 3; 1,2];
fu = zeros(Nu,1);% the right hand side of u
g = zeros(Np,1); % the right hand side of p
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end
if ~isfield(option,'fquadorder')
    option.fquadorder = 4;   % default order is 3
end
if ~isempty(pde.f)
    elem2dofu = [elem2edge NE+elem2edge 2*NE+(1:NT)'];
    % quadrature points in the barycentric coordinate
    [lambda,w] = quadpts(option.fquadorder);
    nQuad = size(lambda,1);
    ft = zeros(NT,7);
    for p = 1:nQuad
        % quadrature points in the x-y coordinate
        pxy = lambda(p,1)*node(elem(:,1),:) ...
            + lambda(p,2)*node(elem(:,2),:) ...
            + lambda(p,3)*node(elem(:,3),:);
        fup = pde.f(pxy);
        for j = 1:3
            j1 = localEdge(j,1); j2 = localEdge(j,2);
            phi_j = lambda(p,j1)*Clambda(:,:,j2)-lambda(p,j2)*Clambda(:,:,j1);
            ft(:,j)   = ft(:,j)  + w(p)*dot(fup,phi_j,2);
            psi_j = lambda(p,j1)*Clambda(:,:,j2)+lambda(p,j2)*Clambda(:,:,j1);
            ft(:,3+j) = ft(:,3+j)+ w(p)*dot(fup,psi_j,2);
        end
        chi = 27*(lambda(p,1)*lambda(p,2)*Clambda(:,:,3) + ...
                  lambda(p,2)*lambda(p,3)*Clambda(:,:,1) + ...
                  lambda(p,3)*lambda(p,1)*Clambda(:,:,2));
        ft(:,7) = ft(:,7) + w(p)*dot(fup,chi,2);
    end
    ft = ft.*repmat(area,1,7);
    fu = accumarray(elem2dofu(:),ft(:),[Nu 1]);
end
clear pxy fup ft phi_j psi_j chi

%% Boundary condition and graddiv part of A
[u,p,ufreeDof,pDof,utbd] = getbdStokesBDM1B;
assembleTime = cputime - t;

%% Solve the system of linear equations
% set up solver type
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if Ndof <= 1e5  % Direct solver for small size systems
        solver = 'direct';
    else             % Multigrid-type  solver for large size systems
        solver = 'mg';
    end
else
    solver = option.solver;
end
% solve the system
t = cputime;
% get submatrices of ufreeDof
A0 = A(ufreeDof,ufreeDof);
B0 = B(:,ufreeDof);
f0 = fu(ufreeDof);
g0 = g;
switch solver
    case 'direct'
    bigA = [A0, B0'; ...
            B0, sparse(Np,Np)];
    bigF = [f0; g0];
    bigu = [u; p];
    bigFreeDof = [ufreeDof; Nu+pDof];
    bigu(bigFreeDof) = bigA(1:end-1,1:end-1)\bigF(1:end-1);
    u = bigu(1:Nu);
    p = bigu(Nu+1:end);
    residual = norm(bigF - bigA*[u(ufreeDof); p]);
    info = struct('solverTime',cputime - t,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'mg'
        option.solver  = 'WCYCLE';
        [u(ufreeDof),p,info] = mgstokesBDM(A0,B0,f0,g0,u,p,node,elemunSort,ufreeDof,option);
    case 'asmg'
        [u(ufreeDof),p,info] = asmgstokes(A0,B0,f0,g0,u,p,node,elemunSort,bdFlag,ufreeDof,option);
end

%% Post-process
if length(pDof)~= Np
    c = sum(p.*area)/sum(area);
    p = p - c;
end
w = R*u + invMvb*utbd;
% drop the bubble part in the output
% u = u(1:2*NE);
w = w(1:N+NE);

%% Output
soln = struct('u',u,'p',p,'w',w);
eqn = struct('A',A0,'B',B0,'Me',Meb,'Mv',Mv,'f',f0,'g',g0,...
             'edge',edge,'ufreeDof',ufreeDof,'pDof',pDof);
info.assembleTime = assembleTime;

%%======================================================================
% subfunction getbdStokesBDM1B
%%======================================================================
    function [u,p,ufreeDof,pDof,utbd] = getbdStokesBDM1B
        %% Boundary condtion of Stokes equations: BDM1B-P0 elements

        % Initial set up
        % Nw = N+NE+NT; Nu = 2*NE+NT; Np = NT; Ndof = Nw+Nu+Np;
        utbd = zeros(Nw,1);
        u = zeros(Nu,1);
        p = zeros(Np,1);
        ufreeDof = (1:Nu)';
        pDof = (1:Np-1)';

        if ~exist('bdFlag','var'), bdFlag = []; end
        if ~isfield(pde,'g_D'), pde.g_D = []; end
        if ~isfield(pde,'g_N'), pde.g_N = []; end
        if ~isfield(pde,'g_R'), pde.g_R = []; end
        if isempty(pde.g_D) && isempty(pde.g_N) && isempty(pde.g_R)
            bdFlag = [];
        end

        % Find Dirichlet boundary dof: fixedDof and pDof
        isFixedEdge = false(NE,1);
        if ~isempty(bdFlag) % if bdFlag is not empty
            % find out the Dirichlet boundary dof
            isDirichlet(elem2edge(bdFlag(:)==1)) = true;
            Dirichlet = edge(isDirichlet,:);
            isFixedEdge(isDirichlet) = true;% dof on D-edges
            fixedEdge = find(isFixedEdge);
            fixedDof = [fixedEdge; fixedEdge+NE];
            ufreeDof = setdiff((1:Nu)',fixedDof);
            % construct the edgeSign
            edgeSign = ones(NE,1);
            idx = (bdFlag(:,1) ~= 0 ) & (elemSign == -1) ; % the first edge is on boundary
            edgeSign(elem2edge(idx,1)) = -1;
            idx = (bdFlag(:,2) ~= 0 ) & (elemSign ==  1) ; % the second edge is on boundary
            edgeSign(elem2edge(idx,2)) = -1;
            idx = (bdFlag(:,3) ~= 0 ) & (elemSign == -1) ; % the third edge is on boundary
            edgeSign(elem2edge(idx,3)) = -1;
        end

        % Compute the boundary integral
        if ~isempty(fixedDof) && ~isempty(pde.g_D) && ~(isnumeric(pde.g_D) && (pde.g_D ==0))
            % 1. Normal component of u is imposed strongly
            if (isnumeric(pde.g_D) && length(pde.g_D) == Nu)
                u(fixedDof) = pde.g_D(fixedDof);
            else
                u(fixedDof) = faceinterpolate(pde.g_D, node,edge(fixedEdge,:),'BDM1');
            end
            % 2. Tangential component of u is imposed weakly
            % 4th order for the line integral
            [lambdagD,wgD] = quadpts1(4);
            nQuadgD = size(lambdagD,1);
            % quadrat = cputime bases 1--3--2
            bdphi = zeros(nQuadgD,3);
            bdphi(:,1) = (2*lambdagD(:,1)-1).*lambdagD(:,1);
            bdphi(:,2) = (2*lambdagD(:,2)-1).*lambdagD(:,2);
            bdphi(:,3) = 4*lambdagD(:,1).*lambdagD(:,2);
            ve = node(Dirichlet(:,2),:) - node(Dirichlet(:,1),:);
            ge = zeros(size(Dirichlet,1),3);
            int_left = zeros(size(Dirichlet,1),2);
            int_right = zeros(size(Dirichlet,1),2);
            int_mid = zeros(size(Dirichlet,1),2);
            for pp = 1:nQuadgD
                ppxy = lambdagD(pp,1)*node(Dirichlet(:,1),:) ...
                     + lambdagD(pp,2)*node(Dirichlet(:,2),:);
                gDp = pde.g_D(ppxy);
                int_left = int_left + wgD(pp)*gDp*bdphi(pp,1);
                int_right = int_right + wgD(pp)*gDp*bdphi(pp,2);
                int_mid = int_mid + wgD(pp)*gDp*bdphi(pp,3); % interior bubble
            end
            ge(:,1) = dot(int_left,ve,2).*edgeSign(fixedEdge);
            ge(:,2) = dot(int_right,ve,2).*edgeSign(fixedEdge);
            ge(:,3) = dot(int_mid,ve,2).*edgeSign(fixedEdge);
            utbd(1:N) = accumarray(Dirichlet(:), [ge(:,1); ge(:,2)],[N,1]);
            utbd(N+fixedEdge) = ge(:,3);
        end

        % Modify the right hand side
        fu = fu- A*u - Meb*(C*(invMvb*utbd));
         g = g - B*u;
         g = g - mean(g);
    end
end

repo: ifem


text:
function item = cfg_branch(varargin)

% This is the branch configuration item class for non-executable
% branches. It implements branch harvest, all_set, get_strings.
%
% Data structure
% ==============
% Description fields
%    * name  - display name of config item
%    * tag   - tag of the menu item
%    * val   - 1xn cell array of cfg_item objects
%    * check - (optional) function handle to implement configuration
%              specific subsasgn checks based on the harvested subtree
%              rooted at this node
%    * help  - help text
% GUI/job manager fields
%    * expanded
%    * hidden
% All fields are inherited from the generic configuration item class.
%
% Public Methods
% ==============
%    * get_strings - returns name of object
%    * gettag      - returns tag
%    * help        - returns help text
%    * harvest     - returns struct, field names correspond to tags of
%                    items in .val field
%    * all_set     - returns all(all_set(item.val{...}))
%
% Output in Job Structure (harvest)
% =================================
% The resulting structure is a struct. Its fieldnames correspond to the
% tags of the cfg_items in item.val, the value of each field is the
% harvested data of the corresponding child item.
%
% The layout of the configuration tree and the types of configuration items
% have been kept compatible to a configuration system and job manager
% implementation in SPM5 (Statistical Parametric Mapping, Copyright (C)
% 2005 Wellcome Department of Imaging Neuroscience). This code has been
% completely rewritten based on an object oriented model of the
% configuration tree.
%
%   The resulting data structure is a struct, with fieldnames according
%   to the 'tag's of the child nodes.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_branch.m 5678 2013-10-11 14:58:04Z volkmar $

rev = '$Rev: 5678 $'; %#ok

myclass = mfilename;
% Get local fields and defaults from private/mysubs_fields
[fn, defs] = mysubs_fields;

if nargin == 1
    if isstruct(varargin{1})
        % assume input is a struct to be converted back into a class
        % return with error if this does not work
        if numel(fieldnames(varargin{1})) == numel(fn)+2 && ...
                all(isfield(varargin{1}, [fn(:)', {'cfg_item' 'cfg_intree'}]))
            gitem = varargin{1}.cfg_item;
            sitem = rmfield(varargin{1},{'cfg_item','cfg_intree'});
            item  = class(sitem, myclass, gitem, cfg_intree);
            return;
        else
            cfg_message('matlabbatch:constructor:reclassify', ['Don''t know how to convert this ' ...
                                'into class ''%s''.'], myclass);
        end;
    end;
    if isa(varargin{1},myclass)
        item = varargin{1};
        return;
    end;
end;

mxpnargin = 4; % Max 4 arguments to parent initialisation
pnargin = min([nargin,mxpnargin]);
switch nargin
    case 0
        gitem = cfg_item;
    case {1,2,3,4,5}
        gitem = cfg_item(varargin{1:pnargin});
    otherwise
        cfg_message('matlabbatch:constructor:nargin', 'Wrong number of arguments.');
end;
item = class(struct([]), myclass, gitem, cfg_intree);
if nargin > mxpnargin
    item.cfg_item.val{1} = varargin{mxpnargin+1};
    mxpnargin = mxpnargin+1;
end;


text:
function [alpha, E, A] = inexact_alm_rlr(X, K, D, tol, maxIter)

% This matlab code implements the inexact augmented Lagrange multiplier
% method for Robust linear regresion with known alpha
%
% alpha- k x n matrix of known coefficients
%
% X - m x n matrix of observations/data (required input)
%
% tol - tolerance for stopping criterion.
%
% maxIter - maximum number of iterations
%
%
% min_D \|E\|_1, s.t. X = D*alpha +E;
%
% The laplacian function: L(A,E,Y,u) =  |E|_1  + \mu/2* |X- D*alpha-E+Y/mu|_F^2;
%
%
% Xianbiao Shu (xshu2@illinois.edu)
% Copyright: Mitsubishi Electric Research Lab
% Reference: X. Shu, F. Porikli, N. Ahujia "Robust Orthonormal Subspace Learning: Efficient Recovery of Corrupted Low-rank Matrices". CVPR 2014;

addpath PROPACK;

[m, n] = size(X);

% initialize
% alpha   % Low-rank coeffs

E = zeros(size(X)); % sparse innovation
A = X;

Z = 0; % Z = Y/mu

norm_inf = norm(X(:), inf);
d_norm = norm(X, 'fro');

mu = 10*5e-2/norm_inf; % this parameter can be tuned
rho = 1.5;       %1.3    % this parameter can be tuned
mu_bar = mu * 1e7;

iter = 0;
Iteration = 0;
converged = false;

while ~converged

    iter = iter + 1;

    E_temp = X - A + Z;

    %% Intensity threshold
    E = max(abs(E_temp) - 1/mu, 0).*sign(E_temp);

    %% Sparse Representation
    A_ini = X -E + Z;

%%   given D and A

    [U, S, V] = svd(D, 0);
    diag_S = diag(S);
    SV = length(diag_S(diag_S>0));

    alpha = V(:, 1:SV)*diag(1./diag_S(1:SV))*(U(:, 1:SV))'*A_ini;
    A = D*alpha;

    Iteration = Iteration + 1;

    Error = X - A - E;
    Z = (Z + Error)/rho;

    mu = min(mu*rho, mu_bar);

    %% stop Criterion
    stopCriterion = norm(Error, 'fro') / d_norm;
    if stopCriterion < tol || iter >= maxIter
        converged = true;
    end

     disp(['#Iter ' num2str(Iteration)  ...
            ' |E|_0 ' num2str(length(find(abs(E)>0)))...
            ' stopCriterion ' num2str(stopCriterion)]);

end