text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [vertices, sph_verts, faces, fvec, new_name]=align_CPS_SHREC(filename, confs)
global fact;
global sgm;

gran = confs.GroupAlpha;  % gran = # of alphas to process together
res = [confs.BaseRes confs.HierarchyStep confs.HierarchyDepth confs.Top_K confs.GammaRes];
% base res R + step of hierarchy Hs + depth of hierarchy Hd + top N + 3rd Angle res (gammares)

switch upper(deblank(char(confs.NormalizeSize)))
    case 'YES'
        scale = 1;
    case 'NO'
        scale = 0;
end

% factorial(170) = Inf
for i=0:170
    fact(i+1) = factorial(i);
end

utl_sgm(15);

% Load a template object
load(confs.Template);
[fvec, max_d] = fixed_fvec(fvec,confs.MaxSPHARMDegree,scale);
%tvertices = vertices; tsph_verts=sph_verts; tfvec = fvec;
tfvec = fvec;

load(filename);
[fvec, max_d] = fixed_fvec(fvec,confs.MaxSPHARMDegree,scale);
[path,name,ext] = fileparts(filename);

if ~exist('faces', 'var') | ~exist('vertices', 'var') | ~exist('sph_verts', 'var') | ~exist('fvec', 'var')
    disp('One or more of faces, vertices, spherical vertices, or SPHARM descriptor are missing');
    return;
end

rmsd_org = SPHARM_rmsd(fvec, tfvec);

if rmsd_org > 0

    % create samples in rotation space
    [alpha, beta, gamma] = utl_eas(res); % euler_angle hierarchy

    % initial alignment using ICP
    [fvec_icp, Robj_icp] = match_icp(fvec, tfvec, max_d);
    [fvec_icp, Aprm] = match_param_hie(fvec_icp,tfvec,alpha,beta,gamma,gran,res,max_d);
    [fvec_icp, Robj_icp] = match_cps(fvec_icp, tfvec, max_d);

    % assume that fvec and atlas are roughly aligned in the object space
    [fvec_cps, Robj_cps] = match_cps(fvec, tfvec, max_d);

    rmsd_icp = SPHARM_rmsd(fvec_icp, tfvec);
    rmsd_cps = SPHARM_rmsd(fvec_cps, tfvec);

    disp(sprintf('RMSD: org %0.3f, icp %0.3f, cps %0.3f',rmsd_org,rmsd_icp,rmsd_cps))

    if rmsd_org > min(rmsd_icp,rmsd_cps)
        if rmsd_icp < rmsd_cps
            fvec = fvec_icp;
            Robj = Robj_icp;
        else
            fvec = fvec_cps;
            Robj = Robj_cps;
        end
    end

    % assume that fvec and atlas are roughly aligned in both object and
    % parameter spaces
    for k = 1:5
        rmsd(1) = SPHARM_rmsd(fvec, tfvec);

        % use cps to align objects together first
        [fvec_a, Robj] = match_cps(fvec, tfvec, max_d);
        rmsd(2) = SPHARM_rmsd(fvec_a, tfvec);

        if rmsd(2)<rmsd(1)
            fvec = fvec_a; % update
        end

        % assume aligned in object space, register the parameterization
        [fvec, Aprm] = match_param_hie(fvec,tfvec,alpha,beta,gamma,gran,res,max_d);
        rmsd(3) = SPHARM_rmsd(fvec, tfvec);

        disp(sprintf('Base Res %d, Step %d, Depth %d, Top %d, Gamma %d: RMSD %0.3f => %0.3f => %0.3f',res,rmsd));

        if sum(abs(Aprm(:)))==0
            break;
        end
    end
else
    disp(sprintf('RMSD = %d: Individual (%s) is the same as template (%s)',rmsd_org,filename,confs.Template));
end

new_name = sprintf('%s/%sSHREC_reg.mat',confs.OutDirectory, name(1:end-3));
if exist(new_name,'file')
    prompt = {'Enter new filename:'};
    dlg_title = 'New File Name';
    num_lines = 1;
    def = {new_name};
    answer = inputdlg(prompt,dlg_title,num_lines,def);
    new_name = answer{1};
end
save(new_name, 'fvec');

clear('fact','sgm');

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [c, ceq] = ma_truConstraint(stateLog, eventID, lbTru, ubTru, bodyIDApply, celBodyData, maData)
%ma_semiMajorAxisConstraint Summary of this function goes here
%   Detailed explanation goes here
    normFact = pi;

    if(ischar(eventID) && strcmpi(eventID,'final'))
        eventNum = max(stateLog(:,13));
    else
%         hMAMainGUI = findall(0,'tag','ma_MainGUI');
%         maData = getappdata(hMAMainGUI,'ma_data');
        [~, eventNum] = getEventByID(eventID, maData.script);
    end

    eventLog = stateLog(stateLog(:,13)==eventNum,:);

    bodyEventLog = eventLog(eventLog(:,8)==bodyIDApply,:);
    if(isempty(bodyEventLog))
        finalEntry = eventLog(end,:);
    else
        finalEntry = bodyEventLog(end,:);
    end

    bodyID = finalEntry(8);

    if(bodyID == bodyIDApply || bodyIDApply==-1)
        bodyInfo = getBodyInfoByNumber(bodyID, celBodyData);
        gmu = bodyInfo.gm;
        rVect = finalEntry(2:4)';
        vVect = finalEntry(5:7)';

        [~, ~, ~, ~, ~, tru] = getKeplerFromState(rVect,vVect,gmu);

        if(lbTru <= 0 && ubTru <= 0)
            lbTru = AngleZero2Pi(lbTru);
            ubTru = AngleZero2Pi(ubTru);
        elseif(lbTru<0 && ubTru >=0)
            if(tru >= AngleZero2Pi(lbTru) && tru <= 2*pi)
                lbTru = AngleZero2Pi(lbTru);
                ubTru = 2*pi;
            elseif(tru >= 0 && tru <= ubTru)
                lbTru = 0;
            end
        end

        if(lbTru == ubTru)
            c = [0 0];
            ceq(1) = tru - ubTru;
        else
            c(1) = lbTru - tru;
            c(2) = tru - ubTru;
            ceq = [0];
        end
        c = c/normFact;
        ceq = ceq/normFact;
    else
        c = [0 0];
        ceq = [0];
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function disp(t,name)
%DISP Command window display of a matricized tensor (tenmat).
%
%   DISP(T) displays a tensor as matrix with no name.
%
%   DISP(T,NAME) display a tensor as matrix with the given name.
%
%   See also TENMAT, TENMAT/DISPLAY.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

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

if ~exist('name','var')
    name = 'ans';
end

fprintf('%s is a matrix corresponding to a tensor of size %s\n',...
        name,tt_size2str(t.tsize));
fprintf('\t%s.rindices = %s (modes of tensor corresponding to rows)\n',...
        name,['[ ' num2str(t.rindices) ' ]']);
fprintf('\t%s.cindices = %s (modes of tensor corresponding to columns)\n',...
        name,['[ ' num2str(t.cindices) ' ]']);

if isempty(t.data)
    fprintf('\t%s.data = []\n',name);
else
    fprintf('\t%s.data = \n',name);
    output = tt_matrix2cellstr(t.data);
    fprintf('\t\t%s\n',output{:});
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% Fig. 6.12   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%

clear all;
%close all;
clf

num1=[10 10];
den1=[1 10];
w=logspace(-2,3,100);
[m1,p1]=bode(num1,den1,w);
num2=[10 -10];
[m2,p2]=bode(num2,den1,w);
figure(1)
loglog(w,m1,w,m2);
xlabel('\omega (rad/sec)');
ylabel('Magnitude');
title('Fig. 6.12 Bode Plot for a NMP System (a) magnitude');
bodegrid;
%pause;
figure(2)
semilogx(w,p1,w,p2,'--');
xlabel('\omega (rad/sec)');
ylabel('Phase (deg)');
title('Fig 6.12 (b) phase');
bodegrid;

author: Sable

repo: mcbench-benchmarks


text:
function [ftps2] = mmps22ftps2(mmps2)
% Convert acceleration from millimeters per square-second to feet per
% second squared
% Chad A. Greene 2012
ftps2 = mmps2*0.00328084;


text:
function [m, o, lambda, eigv]  = mean(o,varargin)
% mean of a list of orientations, principle axes and moments of inertia
%
% Syntax
%   [m, q, lambda, V] = mean(ori)
%   [m, q, lambda, V] = mean(ori,'robust')
%   [m, q, lambda, V] = mean(ori,'weights',weights)
%
% Input
%  ori      - list of @orientation
%
% Options
%  weights  - list of weights
%
% Output
%  m      - mean @orientation
%  o      - crystallographic equivalent @orientation projected to fundamental region
%  lambda - principle moments of inertia
%  V      - principle axes of inertia (@orientation)
%
% See also
% BinghamODF

if isempty(o)
  m = o;
  m.a = NaN; m.b = NaN; m.c = NaN; m.d = NaN; m.i = false;
  if nargout > 2, lambda = zeros(1,4); end
  if nargout > 3, eigv = eye(4); end
  return
elseif length(o) == 1
  m = o;
  if nargout > 1
    eigv = eye(4);
    lambda = [1,0,0,0];
  end
  return;
end

if check_option(varargin,'noSymmetry')

  [m, lambda, eigv] = mean@quaternion(o,varargin{:});

else

  s = size(o);

  % first approximation
  m = get_option(varargin,'q0');
  if isempty(m), m = o.subSet(find(~isnan(o.a),1)); end

  % project around q_mean
  o = project2FundamentalRegion(o,m);

  % compute mean without symmetry
  [m, lambda, eigv] = mean@quaternion(o,varargin{:});

  d = abs(quat_dot(o,m));
  if min(d(:)) < cos(10*degree)
    o = project2FundamentalRegion(o,m);
    [m, lambda, eigv] = mean@quaternion(o,varargin{:});
  end

  if nargout > 1, o = reshape(project2FundamentalRegion(o,m),s); end

end
m.i = false;

author: mtex-toolbox

repo: mtex

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


text:
%=========================================================
%
%This is a prog in the MutualInfo 0.9 package written by
% Hanchuan Peng.
%
%Disclaimer: The author of program is Hanchuan Peng
%      at <penghanchuan@yahoo.com> and <phc@cbmv.jhu.edu>.
%
%The CopyRight is reserved by the author.
%
%Last modification: April/19/2002
%
%========================================================
%
%
%batch compiling of the mex files from cpp src
%by Hanchuan Peng
%April/16/2002

list = dir('*.cpp');

for i=1:length(list),

  fprintf('building mex(dll) of %s\n',list(i).name);
  mex(list(i).name);

end;

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function [tar,non] = get_tar_non(scr,key)
% Divides scores into target and non-target scores using
% information in a key.
% Inputs:
%   scr: A Scores object.
%   key: A Key object.
% Outputs:
%   tar: A vector of target scores.
%   non: A vector of non-target scores.

assert(nargin==2)
assert(isa(scr,'Scores'))
assert(isa(key,'Key'))
assert(scr.validate())
assert(key.validate())

scr = scr.align_with_ndx(key);
scores = scr.scoremat;
tarndx = key.tar & scr.scoremask;
nonndx = key.non & scr.scoremask;
tar = scores(tarndx(:))';
non = scores(nonndx(:))';
end

author: nesl


text:
function hole_point = p07_hole_point ( hole_index, m )

%*****************************************************************************80
%
%% P07_HOLE_POINT returns a point inside a given hole in problem 7.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer HOLE_INDEX, the index of the hole.
%
%    Input, integer M, the spatial dimension.
%
%    Output, real HOLE_POINT(M), a point in the hole
%
  hole_point(1:2) = [ ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function a = r8vec_to_r8cb ( m, n, ml, mu, x )

%*****************************************************************************80
%
%% R8VEC_TO_R8CB copies an R8VEC into a R8CB matrix.
%
%  Discussion:
%
%    In C++ and FORTRAN, this routine is not really needed.  In MATLAB,
%    a data item carries its dimensionality implicitly, and so cannot be
%    regarded sometimes as a vector and sometimes as an array.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns in the array.
%
%    Input, integer ML, MU, the lower and upper bandwidths.
%
%    Input, real X((ML+MU+1)*N), the vector to be copied into the array.
%
%    Output, real A(ML+MU+1,N), the array.
%
  for j = 1 : n
    for i = 1 : ml + mu + 1

      if ( 1 <= i + j - mu - 1 & i + j - mu - 1 <= m )
        a(i,j) = x(i+(ml+mu+1)*(j-1));
      else
        a(i,j) = 0.0;
      end

    end
  end

  return
end


text:
% run_simSteamShaped_fast.m
% Jamie Near, Sunnybrook Research Institute 2022.
%
% USAGE:
% out=run_simSteamShaped_fast(spinSys);
%
% DESCRIPTION:
% This script simulates a STEAM experiment with fully shaped RF
% pulses. Coherence order filtering is employed to eliminate unwanted
% coherences.  Furthermore, simulations are run at various locations in
% space to account for the within-voxel spatial variation of the metabolite
% signal.  Summation across spatial positions is performed. The MATLAB
% parallel computing toolbox (parfor loop) was used to accelerate the
% simulations.  Acceleration is currently performed in the direction of the
% slice selective pulse along the x-direction, but this can be changed.  Up
% to a factor of 12 acceleration can be achieved using this approach.  If
% the parallel processing toolbox is not available, then replace
% the "parfor" loop with a "for" loop.
%
% INPUTS:
% To run this script, there is technically only one input argument:
% spinSys           = spin system to simulate
%
% However, the user should also edit the following parameters as
% desired before running the function:
% RFWaveform        = name of rf pulse waveform used for both (2nd & 3rd) selective 90 degree pulses.
% Tp                = duration of rf pulses[ms]
% Bfield            = Magnetic field strength in [T]
% Npts              = number of spectral points
% sw                = spectral width [Hz]
% Bfield            = magnetic field strength [Tesla]
% lw                = linewidth of the output spectrum [Hz]
% thkX              = slice thickness of x refocusing pulse [cm]
% thkY              = slice thickness of y refocusing pulse [cm]
% fovX              = full simulation FOV in the x direction [cm]
% fovY              = full simulation FOV in the y direction [cm]
% nX                = number of spatial grid points to simulate in x-direction
% nY                = number of spatial grid points to simulate in y-direction
% taus              = vector of pulse sequence timings  [ms]
%
% OUTPUTS:
% out               = Simulation results, summed over all space.

function out=run_simSteamShaped_fast(sys)
tic
% ************INPUT PARAMETERS**********************************
rfWaveform='ex40.b4_384_14.pta'; %name of RF pulse waveform.
Tp=1.920; %duration of RF pulses[ms]
flipAngle=90;  %Flip Angle of the RF pulses [degrees]
Npts=8192; %number of spectral points
sw=6000; %spectral width [Hz]
Bfield=2.89; %magnetic field strength [Tesla]
lw=1; %linewidth of the output spectrum [Hz]
thkX=2.5; %slice thickness of x RF pulse [cm]
thkY=2.5; %slice thickness of y RF pulse [cm]
fovX=5; %size of the full simulation Field of View in the x-direction [cm]
fovY=5; %size of the full simulation Field of View in the y-direction [cm]
nX=48; %Number of grid points to simulate in the x-direction
nY=48; %Number of grid points to simulate in the y-direction
tau1=6; %TE for STEAM sequence [ms]
tau2=32; %TM for STEAM sequence [ms]
centreFreq=2.3; %Centre frequency of simulation [ppm]
% ************END OF INPUT PARAMETERS**********************************

%set up spatial grid
x=linspace(-fovX/2,fovX/2,nX); %X positions to simulate [cm]
y=linspace(-fovY/2,fovY/2,nY); %y positions to simulate [cm]

%Load RF waveform
RF=io_loadRFwaveform(rfWaveform,'exc',0);

gamma=42577000; %gyromagnetic ratio

%Load spin systems
%load spinSystems
%sys=eval(['sys' spinSys]);

%If length of RF pulse is >200 pts, resample to 100 pts to reduce
%computational workload
if size(RF.waveform,1)>200
    RF=rf_resample(RF,100);
end

Gx=(RF.tbw/(Tp/1000))/(gamma*thkX/10000); %[G/cm]
Gy=(RF.tbw/(Tp/1000))/(gamma*thkY/10000); %[G/cm]

%Initialize structures:
d_temp=cell(1,1);
d=cell(1,1);

%loop through space: If you are using the parfor loops below, and you are
%using an older version of MATLAB (e.g.R2012), don't forget to initialize
%the parallel processing toolbox workers using 'matlabpool open N' (for N
%workers, 12 max).  I don't think this is necessary for newer versions of
%MATLAB.

%First loop through all x-positions, simulating only the first refocusing
%pulse.
%First loop through x direction (first refoc pulse only);

%for X=1:length(x)  %Use this if you don't have the MATLAB parallel processing toolbox
parfor X=1:length(x)  %Use this if you have the MATLAB parallel processing toolbox
        disp(['Executing X-position ' num2str(X) ' of ' num2str(length(x)) '!!!']);
        d_temp{X}=sim_steam_shaped_fastRF1(Bfield,sys,tau1,tau2,RF,Tp,x(X),Gx,flipAngle,centreFreq);
end

%calculate the average density matrix (Doing this inside a separate for
%loop because I couldn't figure out how to do this inside the parfor loop):
for X=1:length(x)
        d{1}=sim_dAdd(d{1},d_temp{X});
end

% %Initialize structures:
out_temp=cell(length(y),1);
out=struct([]);

%Now loop through y direction (second refoc pulse only);
%for Y=1:length(y) %Use this if you don't have the MATLAB parallel processing toolbox
parfor Y=1:length(y) %Use this if you do have the MATLAB parallel processing toolbox
%            disp(['Executing Y-position ' num2str(Y) ' of ' num2str(length(y)) '!!!']);
            out_temp{Y}=sim_steam_shaped_fastRF2(d{1},Npts,sw,Bfield,lw,sys,tau1,tau2,...
                RF,Tp,y(Y),Gy,flipAngle,centreFreq);
end

%Now combine the outputs;  Again, doing this inside a separate for loop
%becuase I can't figure out how to do this inside the parfor loop:
for Y=1:length(y)
            out=op_addScans(out,out_temp{Y});
end

%For consistent scaling across different shaped simulations, we need to :
%1.  Scale down by the total number of simulations run (since these were
%    all added together.
numSims=(nX*nY);
out=op_ampScale(out,1/numSims);

%2.  Scale by the total size of the simulated region, relative to the size
%    of the voxel.
voxRatio=(thkX*thkY)/(fovX*fovY);
out=op_ampScale(out,1/voxRatio);

toc
end

%Nested Function #1
function d = sim_steam_shaped_fastRF1(Bfield,sys,TE,TM,RF,tp,dx,Gx,flipAngle,centreFreq)
%
% USAGE:
% d = sim_steam_shaped_fastRF1(Bfield,sys,TE,TM,RF,tp,dx,Gx,flipAngle,centreFreq)
%
% DESCRIPTION:
% This function simulates only the first bit of the STEAM experiment, up to
% the beginning of the third RF pulse.  The excitation is
% simulated as an instantaneous rotation, and the refocusing pulse is
% simulated as a shaped rotation.
%
% This code is designed to be used in highly-accelerated shaped simulations,
% using the method described by Yan Zhang et al. Med Phys 2017;44(8):
% 4169-78.
%
% This code employs coherence order filtering to select only desired
% signal.  Finally, this code simulates the spectrum at a given point in
% space (x), given the values of the slice selection gradient (Gx).  In
% order to fully simulate the STEAM experiment, you have to run this
% simulation many times at various points in space (x), followed by
% sim_steam_shaped_fastRF2.m, at all points in space (y).
%
% INPUTS:
% Bfield    = main magnetic field strength in [T]
% sys       = spin system definition structure
% TE        = echo time in [ms].
% TM        = mixing time in [ms].
% RF        = RF pulse definition structure for RF pulses (obtain using 'io_loadRFwaveform.m')
% tp        = RF pulse duration in [ms]
% dx        = position offset in x-direction (corresponding to 2nd STEAM pulse) [cm]
% Gx        = gradient strength for 2nd STEAM pulse [G/cm]
% flipAngle = flip angle of refocusing pulses [degrees] (Optional.  Default = 90 deg)
%
% OUTPUTS:
% out       = simulated spectrum, in FID-A structure format, using STEAM
%             sequence.
if nargin<10
    centreFreq=2.3;
    if nargin<9
        flipAngle=90;
    end
end

%In the steam sequence, it can be common to use an asymmetric RF pulse for the
%90-degree pulse waveform.  In this case, it is conventional for the 2nd
%and 3rd rf pulses to be time-reversed versions of eachother, with the 2nd
%pulse (RF1) being a max-phase pulse, and the 3rd pulse (RF2) being a
%min-phase pulse (i.e. the long tails of both pulse occurring during the TM
%period so that the TE is minimized).  Here, check if the RF pulse is
%asymmetric and if so, make sure that the 2nd and 3rd pulses are max-phase
%and min-phase, respectively:
if RF.rfCentre>0.5
    RF1=rf_timeReverse(RF);
    RF2=RF;
else
    RF1=RF;
    RF2=rf_timeReverse(RF);
end

%Check that the TE and TM values are not too short
if TE<(RF1.rfCentre*tp*2)
    error('ERROR:  TE cannot be less than duration of RF pulse! ABORTING!!');
end
if TM<(RF2.rfCentre*tp*2)
    error('ERROR:  TM cannot be less than duration of RF pulse! ABORTING!!');
end

%Set centre frequency
for k=1:length(sys)
    sys(k).shifts=sys(k).shifts-centreFreq;
end

%Calculate Hamiltonian matrices and starting density matrix.
[H,d]=sim_Hamiltonian(sys,Bfield);

%Calculate new delays by subtracting the pulse duration from tau1 and tau2;
delays=zeros(2);
delays(1)=TE-(RF1.rfCentre*tp*2);
delays(2)=TM-(RF2.rfCentre*tp*2);
if sum(delays<0)
    error(['ERROR! The following taus are too short: ' num2str(find(delays<0)) '.']);
end

%BEGIN PULSE SEQUENCE************
d=sim_excite(d,H,'x');                                    %EXCITE
d=sim_COF(H,d,1);                                       %Keep only +1-order coherences
d=sim_evolve(d,H,delays(1)/2000);                         %Evolve by delays(1)/2
d=sim_gradSpoil(d,H,[Gx,0,0],[dx,0,0],tp*RF1.rfCentre);     %Prewind gradient for 2nd 90 degree pulse (Not sure why, but this only works when Gx is positive.  Intiutively, Gx amplitude should be the opposite of the slice select gradient (i.e. -Gx), but this does not seem to work).
d=sim_shapedRF(d,H,RF1,tp,flipAngle,90,dx,Gx);             %1st shaped 90 degree selective pulse
d=sim_COF(H,d,0);                                       %Keep only 0-order coherences
d=sim_evolve(d,H,(delays(2))/1000);                       %Evolve by delays(2)
%END PULSE SEQUENCE**************

%After running this many times along x, the density matrices should be
%averaged, and then the average density matrix should be passed through
%'sim_press_shaped_fastRef2' at various different y-positions.

end

%Nested Function #2
function out = sim_steam_shaped_fastRF2(d,n,sw,Bfield,linewidth,sys,TE,TM,RF,tp,dy,Gy,flipAngle,centreFreq)
%
% USAGE:
% out = sim_steam_shaped_fastRF(d,n,sw,Bfield,linewidth,sys,TE,TM,RF,tp,dy,Gy,flipAngle,centreFreq)
%
% DESCRIPTION:
% This function simulates only the last bit of the STEAM experiment, from the
% the beginning of the third STEAM pulse, to the end.  The RF
% pulse is simulated as a shaped rotation.
%
% This code is designed to be used in highly-accelerated shaped simulations,
% using the method described by Yan Zhang et al. Med Phys 2017;44(8):
% 4169-78.

% This code employs coherence order filtering to select only desired
% signal.
%
% Finally, this code simulates the spectrum at a given point in space (y),
% given the values of the slice selection gradient (Gy).  In order
% to fully simulate the STEAM experiment, you have to first run
% sim_steam_shaped_fastRF1.m at all points in space (x), followed by
% this code, at all points in space (y).
%
% INPUTS:
% d         = starting density matrix (obtained using 'sim_steam_shaped_fastRF1.m')
% n         = number of points in fid/spectrum
% sw        = desired spectral width in [Hz]
% Bfield    = main magnetic field strength in [T]
% linewidth = linewidth in [Hz]
% sys       = spin system definition structure
% TE        = echo time in [ms].
% TM        = mixing time in [ms].
% RF        = RF pulse definition structure for selective RF pulses (obtain using 'io_loadRFwaveform.m')
% tp        = RF pulse duration in [ms]
% dx        = position offset in x-direction (corresponding to first refocusing pulse) [cm]
% dy        = position offset in y-direction (corresponding to second refocusing pulse) [cm]
% Gx        = gradient strength for 2nd STEAM RF pulse [G/cm]
% Gy        = gradient strength for 3rd STEAM RF pulse [G/cm]
% flipAngle = flip angle of RF pulses [degrees] (Optional.  Default = 90 deg)
%
% OUTPUTS:
% out       = simulated spectrum, in FID-A structure format, using STEAM
%             sequence.

if nargin<14
    centreFreq=2.3;
    if nargin<13
        flipAngle=90;
    end
end

%In the steam sequence, it can be common to use an asymmetric RF pulse for the
%90-degree pulse waveform.  In this case, it is conventional for the 2nd
%and 3rd rf pulses to be time-reversed versions of eachother, with the 2nd
%pulse (RF1) being a max-phase pulse, and the 3rd pulse (RF2) being a
%min-phase pulse (i.e. the long tails of both pulse occurring during the TM
%period so that the TE is minimized).  Here, use the asymmetry factor of
%the pulse to make sure that the 2nd and 3rd pulses are max-phase and
%min-phase, respectively:
if RF.rfCentre>0.5
    RF1=rf_timeReverse(RF);
    RF2=RF;
else
    RF1=RF;
    RF2=rf_timeReverse(RF);
end

%Check that the TE and TM values are not too short
if TE<(RF1.rfCentre*tp*2)
    error('ERROR:  TE cannot be less than duration of RF pulse! ABORTING!!');
end
if TM<(RF2.rfCentre*tp*2)
    error('ERROR:  TM cannot be less than duration of RF pulse! ABORTING!!');
end

%Set centre frequency
for k=1:length(sys)
    sys(k).shifts=sys(k).shifts-centreFreq;
end

%Calculate Hamiltonian matrices and starting density matrix.
[H]=sim_Hamiltonian(sys,Bfield);

%Calculate new delays by subtracting the pulse duration from TE and TM;
delays=zeros(2);
delays(1)=TE-(RF1.rfCentre*tp*2);
delays(2)=TM-(RF2.rfCentre*tp*2);
if sum(delays<0)
    error(['ERROR! The following delays are too short: ' num2str(find(delays<0)) '.']);
end

%BEGIN PULSE SEQUENCE************
d=sim_shapedRF(d,H,RF2,tp,flipAngle,90,dy,Gy);          %2nd shaped 90 degree selective pulse
d=sim_gradSpoil(d,H,[0,Gy,0],[0,dy,0],tp*RF1.rfCentre);   %Rewind gradient for 3rd 90 degree pulse (Not sure why, but this only works when Gy is positive.  Intiutively, Gx amplitude should be the opposite of the slice select gradient (i.e. -Gy), but this does not seem to work).
d=sim_COF(H,d,-1);                                      %Keep only -1 coherences
d=sim_evolve(d,H,delays(1)/2000);                            %Evolve by delays(1)/2
[out,~]=sim_readout(d,H,n,sw,linewidth,90);      %Readout along y (90 degree phase);
%END PULSE SEQUENCE**************

%Correct the ppm scale:
out.ppm=out.ppm-(4.65-centreFreq);

%Fill in structure header fields:
out.seq='steam';
out.te=TE;
out.tm=TM;
out.sim='shaped';

%Additional fields for compatibility with FID-A processing tools.
out.sz=size(out.specs);
out.date=date;
out.dims.t=1;
out.dims.coils=0;
out.dims.averages=0;
out.dims.subSpecs=0;
out.dims.extras=0;
out.averages=1;
out.rawAverages=1;
out.subspecs=1;
out.rawSubspecs=1;
out.flags.writtentostruct=1;
out.flags.gotparams=1;
out.flags.leftshifted=0;
out.flags.filtered=0;
out.flags.zeropadded=0;
out.flags.freqcorrected=0;
out.flags.phasecorrected=0;
out.flags.averaged=1;
out.flags.addedrcvrs=1;
out.flags.subtracted=1;
out.flags.writtentotext=0;
out.flags.downsampled=0;
out.flags.isFourSteps=0;

end

author: CIC-methods

repo: FID-A

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


text:
classdef EigModes < handle

    properties (Access = public)

    end

    properties (Access = private)
        mesh
        dim
        V
        v1
        v2
        D
        eigModesPlotter
        lambda
        freeNodes
        mode1Disp
        mode2Disp
    end

    properties (Access = private)
        designVariable
        stiffnessMatComputer
        bendingMatComputer
    end

    methods (Access = public)

        function obj = EigModes(cParams)
            obj.init(cParams)
            obj.createEigModesPlotter();
        end

        function plot(obj,A,iter)
            p = obj.eigModesPlotter;
            [m1,m2] = obj.computeBucklingModes(obj.v1,obj.v2);
            obj.mode1Disp = m1;
            obj.mode2Disp = m2;
            p.plot(A,m1,m2,iter,obj.D)
        end

        function fx = provideFunction(obj,eigNum)
            obj.computeEigenModesAndValues();
            obj.lambda = obj.computeLambda();
            gamma = obj.designVariable.getFirstEigenMode();
            fx = gamma-obj.lambda(eigNum);
        end

       function grad = provideDerivative(obj,eigNum)
            obj.reorderModes(obj.lambda,obj.V,obj.D);
            Belem =  obj.bendingMatComputer.elementalBendingMatrix;
            x = obj.designVariable.getColumnArea();
            nElem = obj.mesh.nelem;
            eigV1 = obj.D(1,1);
            eigV2 = obj.D(2,2);
            difEigs = abs(eigV2-eigV1);
            if difEigs > 1
                dfdx = obj.computeSimpleEig(Belem,x);
            else
                dfdx = obj.computeDoubleEig(Belem,x);
            end
            dfdx(1,nElem+1) = 1;
            dfdx(2,nElem+1) = 1;
            grad = dfdx(eigNum,:);
        end

    end

    methods (Access = private)

        function dfdx = computeSimpleEig(obj,Belem,x)
            d = obj.dim;
            free = obj.freeNodes;
            ndofe = d.ndofPerElement;
            ndofn = d.ndimField;
            W = zeros(d.ndof,2);
            W(free,1) = obj.v1;
            W(free,2) = obj.v2;
            nElem = obj.mesh.nelem;

            for i = 1:nElem
                index = ndofn*(i-1)+1: ndofn*(i-1)+ndofe;
                dx = 2*x(i,1);
                dfdx(1,i) = -dx*(W(index,1)'*Belem(:,:,i)*W(index,1));
                dfdx(2,i) = -dx*(W(index,2)'*Belem(:,:,i)*W(index,2));

            end
%             q = Quadrature.set(obj.mesh.type);
%             q.computeQuadrature('LINEAR');
%             l = sum(obj.mesh.computeDvolume(q));
%             dfdx(1,:) = 1./1.*dfdx(1,:);
%             dfdx(2,:) = 1./1.*dfdx(2,:);
        end

        function dfdx = computeDoubleEig(obj,Belem,x)
            d    = obj.dim;
            free = obj.freeNodes;
            ndofe = d.ndofPerElement;
            ndofn = d.ndimField;
            nElem = obj.mesh.nelem;
            W1    = zeros(d.ndof,1);
            W2    = zeros(d.ndof,1);
            dW1   = zeros(nElem,1);
            dW2   = zeros(nElem,1);
            dW1W2 = zeros(nElem,1);
            W1(free,1) = obj.v1;
            W2(free,1) = obj.v2;
            for i=1:nElem
                index = ndofn*(i-1)+1: ndofn*(i-1)+ndofe;
                dx = 2*x(i,1);
                dW1(i,1)= dx*(W1(index,1)'*Belem(:,:,i)*W1(index,1));
                dW2(i,1)= dx*(W2(index,1)'*Belem(:,:,i)*W2(index,1));
                dW1W2(i,1)= (2*x(i,1))*(W1(index,1)'*Belem(:,:,i)*W2(index,1));
                A = [dW1(i,1) dW1W2(i,1); dW1W2(i,1) dW2(i,1)];
                [U,R] = eigs(A,2,'SM');
                S = sort(diag(R));
                dfdx(1,i) = -S(1);
                dfdx(2,i) = -S(2);
            end
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.mesh                 = cParams.mesh;
            obj.dim                  = cParams.dim;
            obj.designVariable       = cParams.designVariable;
            obj.stiffnessMatComputer = cParams.stiffnessMatComputer;
            obj.bendingMatComputer   = cParams.bendingMatComputer;
        end

        function computeEigenModesAndValues(obj)
            obj.stiffnessMatComputer.compute();
            obj.bendingMatComputer.compute();
            [Kfree,free]  = obj.stiffnessMatComputer.provideFreeStiffnessMatrix();
            obj.freeNodes = free;
            Bfree  = obj.bendingMatComputer.provideFreeBendingMatrix();
            obj.computeEigenFunctionAndValues(Bfree,Kfree);
        end

        function createEigModesPlotter(obj)
            s.mesh = obj.mesh;
            p = EigModesPlotter(s);
            obj.eigModesPlotter = p;
        end

        function l = computeLambda(obj)
            l = sort(diag(obj.D));
        end

        function computeEigenFunctionAndValues(obj,B,K)
            [v,d] = eigs(B,K,2,'SM');
            obj.V  = v;
            obj.D  = d;
        end

        function reorderModes(obj,lambda,V,D)
            if lambda(1)==D(1,1)
                V1=V(:,1);
                V2=V(:,2);
            else
                V1=V(:,2);
                V2=V(:,1);
            end
            obj.v1 = V1;
            obj.v2 = V2;
        end

        function [m1disp,m2disp] = computeBucklingModes(obj,v1,v2)
            N = obj.mesh.nelem;
            Mode1=zeros(2*(N+1),1);
            Mode2=zeros(2*(N+1),1);
            for i=3:2*N
                Mode1(i)=v1(i-2);
                Mode2(i)=v2(i-2);
            end
            m1 = Mode1;
            m2 = Mode2;
            m1disp = m1(1:2:end);
            m2disp = m2(1:2:end);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [Output, ChannelFile, FileFormat] = import_channel(iStudies, ChannelFile, FileFormat, ChannelReplace, ChannelAlign, isSave, isFixUnits, isApplyVox2ras, RefMriFile)
% IMPORT_CHANNEL: Imports a channel file (definition of the sensors).
%
% USAGE:  BstChannelFile = import_channel(iStudies=none, ChannelFile=[ask], FileFormat, ChannelReplace=1, ChannelAlign=[ask], isSave=1, isFixUnits=[ask], isApplyVox2ras=[ask], RefMriFile=[])
%
% INPUT:
%    - iStudies       : Indices of the studies where to import the ChannelFile
%    - ChannelFile    : Full filename of the channels list to import (default: asked to the user)
%    - FileFormat     : Format of the input file ChannelFile
%    - ChannelReplace : 0, do not replace if channel file already exist
%                       1, replace old channel file after user confirmation  (default)
%                       2, replace old channel file without user confirmation
%    - ChannelAlign   : 0, do not perform automatic headpoints-based alignment
%                       1, perform automatic alignment after user confirmation  (default)
%                       2, perform automatic alignment without user confirmation
%    - isSave         : If 1, save the new channel file in the target study
%    - isFixUnits     : If 1, tries to convert the distance units to meters automatically
%                       If 0, does not fix the distance units
%                       If [], ask for the scaling to apply
%    - isApplyVox2ras : If 0, does not use the voxel=>subject transformation
%                       If 1, uses the existing voxel=>subject transformation from the MRI file, if available
%                       If 2, uses the existing voxel=>subject transformation AND the coregistration transformation (see process_import_bids)
%                       If [], ask for user decision
%    - RefMriFile     : Relative file name to a MRI file imported for the current subject
%                       When isApplyVox2ras=1: use this file instead of the reference MRI for getting the vox2ras transformation
%                       When isApplyVox2ras=2: uses this file AND reverts the registration matrix that was possibly
%                                              applied after the volume was imported (to matche the original .nii)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2008-2023

%% ===== PARSE INPUTS =====
Output = [];
if (nargin < 9) || isempty(RefMriFile)
    RefMriFile = [];
end
if (nargin < 8) || isempty(isApplyVox2ras)
    isApplyVox2ras = [];
end
if (nargin < 7) || isempty(isFixUnits)
    isFixUnits = [];
end
if (nargin < 6) || isempty(isSave)
    isSave = 1;
end
if (nargin < 5) || isempty(ChannelAlign)
    ChannelAlign = [];
end
if (nargin < 4) || isempty(ChannelReplace)
    ChannelReplace = 1;
end
if (nargin < 3) || isempty(ChannelFile) || isempty(FileFormat)
    ChannelFile = [];
end
if (nargin < 1) || isempty(iStudies)
    iStudies = [];
end

%% ===== SELECT CHANNEL FILE =====
% If file to load was not defined : open a dialog box to select it
if isempty(ChannelFile)
    isInteractive = 1;
    % Get default import directory and formats
    LastUsedDirs = bst_get('LastUsedDirs');
    DefaultFormats = bst_get('DefaultFormats');
    % Get MRI file
    [ChannelFile, FileFormat] = java_getfile('open', ...
            'Import channel file...', ...              % Window title
            LastUsedDirs.ImportChannel, ...        % Last used directory
            'single', 'files_and_dirs', ...        % Selection mode
            bst_get('FileFilters', 'channel'), ... % File filters
            DefaultFormats.ChannelIn);             % Default ASCII XYZ
    % If no file was selected: exit
    if isempty(ChannelFile)
        return
    end
    % Save default import directory
    LastUsedDirs.ImportChannel = bst_fileparts(ChannelFile);
    bst_set('LastUsedDirs', LastUsedDirs);
    % Save default import format
    DefaultFormats.ChannelIn = FileFormat;
    bst_set('DefaultFormats',  DefaultFormats);
else
    isInteractive = 0;
end

%% ===== LOAD CHANNEL FILE =====
ChannelMat = [];
% Progress bar
isProgressBar = bst_progress('isVisible');
if ~isProgressBar
    bst_progress('start', 'Import channel file', ['Loading file "' ChannelFile '"...']);
end
% Get the file extenstion
[fPath, fBase, fExt] = bst_fileparts(ChannelFile);
if ~isempty(fExt)
    fExt = lower(fExt(2:end));
end
% Import options
ImportOptions = db_template('ImportOptions');
ImportOptions.EventsMode = 'ignore';
ImportOptions.DisplayMessages = 0;
% Load file
switch (FileFormat)
    % ===== MEG/EEG =====
    case 'CTF'
        ChannelMat = in_channel_ctf(ChannelFile);
        FileUnits = 'm';
    case {'FIF', '4D', 'KIT', 'BST-BIN', 'KDF', 'RICOH', 'ITAB', 'MEGSCAN-HDF5'}
        [sFile, ChannelMat] = in_fopen(ChannelFile, FileFormat, ImportOptions);
        if isempty(ChannelMat)
            bst_progress('stop');
            return;
        end
        FileUnits = 'm';
    case 'BST'
        ChannelMat = in_bst_channel(ChannelFile);
        FileUnits = 'm';

    % ===== EEG ONLY =====
    case {'BIDS-SCANRAS-MM', 'BIDS-MNI-MM', 'BIDS-ACPC-MM', 'BIDS-ALS-MM', 'BIDS-CAPTRAK-MM'}
        ChannelMat = in_channel_bids(ChannelFile, 0.001);
        FileUnits = 'mm';
    case {'BIDS-SCANRAS-CM', 'BIDS-MNI-CM', 'BIDS-ACPC-CM', 'BIDS-ALS-CM', 'BIDS-CAPTRAK-CM'}
        ChannelMat = in_channel_bids(ChannelFile, 0.01);
        FileUnits = 'cm';
    case {'BIDS-SCANRAS-M', 'BIDS-MNI-M', 'BIDS-ACPC-M', 'BIDS-ALS-M', 'BIDS-CAPTRAK-M'}
        ChannelMat = in_channel_bids(ChannelFile, 1);
        FileUnits = 'm';

    case 'BESA' % (*.sfp;*.elp;*.eps/*.ela)
        switch (fExt)
            case 'sfp'
                ChannelMat = in_channel_ascii(ChannelFile, {'Name','-Y','X','Z'}, 0, .01);
                ChannelMat.Comment = 'BESA channels';
            case 'elp'
                ChannelMat = in_channel_ascii(ChannelFile, {'Name','Y','X','Z'}, 0, .01);
                ChannelMat.Comment = 'BESA channels';
            case {'eps','ela'}
                ChannelMat = in_channel_besa_eps(ChannelFile);
        end
        FileUnits = 'cm';

    case 'BRAINVISION'
        ChannelMat = in_channel_brainvision(ChannelFile);
        FileUnits = 'm';

    case 'CARTOOL' % (*.els;*.xyz)
        switch (fExt)
            case 'els'
                ChannelMat = in_channel_cartool_els(ChannelFile);
                FileUnits = 'mm';
            case 'xyz'
                ChannelMat = in_channel_ascii(ChannelFile, {'-Y','X','Z','Name'}, 1, .001);
                ChannelMat.Comment = 'Cartool channels';
                FileUnits = 'cm';
        end

    case 'CURRY' % (*.res;*.rs3;*.pom)
        switch (fExt)
            case 'res'
                ChannelMat = in_channel_ascii(ChannelFile, {'%d','-Y','X','Z','%d','Name'}, 0, .001);
                ChannelMat.Comment = 'Curry channels';
            case 'rs3'
                ChannelMat = in_channel_curry_rs3(ChannelFile);
            case 'pom'
                ChannelMat = in_channel_curry_pom(ChannelFile);
        end
        FileUnits = 'mm';

    case 'EEGLAB' % (*.ced;*.xyz)
        switch (fExt)
            case 'ced'
                ChannelMat = in_channel_ascii(ChannelFile, {'indice','Name','%f','%f','X','Y','Z','%f','%f','%f'}, 1, .0875); % Convert normalized coord => average head radius
                ChannelMat.Comment = 'EEGLAB channels';
            case 'xyz'
                ChannelMat = in_channel_ascii(ChannelFile, {'indice','-Y','X','Z','Name'}, 0, .01);
                ChannelMat.Comment = 'EEGLAB channels';
            case 'set'
                ChannelMat = in_channel_eeglab_set(ChannelFile, isFixUnits);
        end
        FileUnits = 'cm';

    case 'EETRAK' % (*.elc)
        ChannelMat = in_channel_ascii(ChannelFile, {'X','Y','Z'}, 3, .001);
        ChannelMat.Comment = 'EETRAK channels';
        FileUnits = 'mm';

    case 'EGI'  % (*.sfp)
        ChannelMat = in_channel_ascii(ChannelFile, {'Name','-Y','X','Z'}, 0, .01);
        ChannelMat.Comment = 'EGI channels';
        FileUnits = 'cm';

    case 'EMSE'  % (*.elp)
        ChannelMat = in_channel_emse_elp(ChannelFile);
        FileUnits = 'm';

    case 'FREESURFER-TSV'
        % Read file
        ChannelMat = in_channel_bids(ChannelFile, 0.001);
        FileUnits = 'm';
        % If we know the destination study: convert from FreeSurfer/Surface coordinates to SCS coordinates
        if ~isempty(iStudies)
            % Get the subject for the first study
            sStudy = bst_get('Study', iStudies(1));
            sSubject = bst_get('Subject', sStudy.BrainStormSubject);
            % Get the subject's MRI
            if isempty(sSubject.Anatomy) || isempty(sSubject.Anatomy(1).FileName)
                error('You need to import the FreeSurfer anatomy before.');
            end
            % Load the MRI
            MriFile = file_fullpath(sSubject.Anatomy(1).FileName);
            sMri = in_mri_bst(MriFile);
            if isempty(sMri) || ~isfield(sMri, 'SCS') || ~isfield(sMri.SCS, 'NAS') || isempty(sMri.SCS.NAS)
                error('Missing fiducials.');
            end
            % Convert coordinates: FreeSurfer surface => MRI
            fcnTransf = @(Loc)bst_bsxfun(@plus, Loc', (size(sMri.Cube(:,:,:,1))/2 + [0 1 0]) .* sMri.Voxsize / 1000)';
            AllChannelMats = channel_apply_transf(ChannelMat, fcnTransf, [], 1);
            ChannelMat = AllChannelMats{1};
            % Convert coordinates: MRI => SCS
            fcnTransf = @(Loc)cs_convert(sMri, 'mri', 'scs', Loc')';
            AllChannelMats = channel_apply_transf(ChannelMat, fcnTransf, [], 1);
            ChannelMat = AllChannelMats{1};
        end

    case {'INTRANAT', 'INTRANAT_MNI'}
        switch (fExt)
            case 'pts'
                ChannelMat = in_channel_ascii(ChannelFile, {'name','X','Y','Z'}, 3, .001);
            case 'csv'
                if strcmpi(FileFormat, 'INTRANAT_MNI')
                    ChannelMat = in_channel_tsv(ChannelFile, 'contact', 'MNI', .001);
                else
                    ChannelMat = in_channel_tsv(ChannelFile, 'contact', 'T1pre Scanner Based', .001);
                end
        end
        ChannelMat.Comment = 'Contacts';
        FileUnits = 'mm';
        [ChannelMat.Channel.Type] = deal('SEEG');

    case 'MEGDRAW'
        ChannelMat = in_channel_megdraw(ChannelFile);
        FileUnits = 'cm';

    case 'LOCALITE'
        ChannelMat = in_channel_ascii(ChannelFile, {'%d','name','X','Y','Z'}, 1, .001);
        ChannelMat.Comment = 'Localite channels';
        FileUnits = 'mm';

    case 'MFF'  % (coordinates.xml)
        [tmp, ChannelMat] = in_fopen_mff(ChannelFile, ImportOptions, 1);
        FileUnits = 'mm';

    case 'NEUROSCAN'  % (*.dat;*.tri;*.txt;*.asc)
        switch (fExt)
            case {'dat', 'txt'}
                ChannelMat = in_channel_neuroscan_dat(ChannelFile);
                FileUnits = 'cm';
            case 'tri'
                ChannelMat = in_channel_neuroscan_tri(ChannelFile);
                FileUnits = 'cm';
            case 'asc'
                ChannelMat = in_channel_neuroscan_asc(ChannelFile);
                FileUnits = 'mm';
        end

    case 'POLHEMUS'  % (*.pos;*.elp)
        switch (fExt)
            case 'pos'
                ChannelMat = in_channel_pos(ChannelFile);
                FileUnits = 'cm';
            case {'pol','txt'}
                ChannelMat = in_channel_ascii(ChannelFile, {'name','X','Y','Z'}, 1, .01);
                ChannelMat.Comment = 'Polhemus';
                FileUnits = 'cm';
            case 'elp'
                ChannelMat = in_channel_emse_elp(ChannelFile);
                FileUnits = 'mm';
        end

    case 'SIMNIBS'
        ChannelMat = in_channel_ascii(ChannelFile, {'%s','X','Y','Z','name'}, 0, .001);
        ChannelMat.Comment = '10-10 electrodes';
        FileUnits = 'mm';

    case 'TVB'
        ChannelMat = in_channel_tvb(ChannelFile);
        FileUnits = 'm';

    case 'XENSOR' % ANT Xensor (*.elc)
        ChannelMat = in_channel_ant_xensor(ChannelFile);
        FileUnits = 'mm';

    case {'ASCII_XYZ', 'ASCII_XYZ_MNI', 'ASCII_XYZ_WORLD'}  % (*.*)
        ChannelMat = in_channel_ascii(ChannelFile, {'X','Y','Z'}, 0, .001);
        ChannelMat.Comment = 'Channels';
        FileUnits = 'mm';
    case {'ASCII_NXYZ', 'ASCII_NXYZ_MNI', 'ASCII_NXYZ_WORLD'}  % (*.*)
        ChannelMat = in_channel_ascii(ChannelFile, {'Name','X','Y','Z'}, 0, .001);
        ChannelMat.Comment = 'Channels';
        FileUnits = 'mm';
    case {'ASCII_XYZN', 'ASCII_XYZN_MNI', 'ASCII_XYZN_WORLD'}  % (*.*)
        ChannelMat = in_channel_ascii(ChannelFile, {'X','Y','Z','Name'}, 0, .001);
        ChannelMat.Comment = 'Channels';
        FileUnits = 'mm';
    case 'ASCII_NXY'  % (*.*)
        ChannelMat = in_channel_ascii(ChannelFile, {'Name','X','Y'}, 0, .000875);
        ChannelMat.Comment = 'Channels';
        FileUnits = 'mm';
    case 'ASCII_XY'  % (*.*)
        ChannelMat = in_channel_ascii(ChannelFile, {'X','Y'}, 0, .000875);
        ChannelMat.Comment = 'Channels';
        FileUnits = '';
    case 'ASCII_NTP'  % (*.*)
        ChannelMat = in_channel_ascii(ChannelFile, {'Name','TH','PHI'}, 0, .0875);
        ChannelMat.Comment = 'Channels';
        FileUnits = '';
    case 'ASCII_TP'  % (*.*)
        ChannelMat = in_channel_ascii(ChannelFile, {'TH','PHI'}, 0, .0875);
        ChannelMat.Comment = 'Channels';
        FileUnits = '';
    otherwise
        error(['File format is not supported: ' FileFormat]);
end
% No data imported
isHeadPoints = isfield(ChannelMat, 'HeadPoints') && ~isempty(ChannelMat.HeadPoints.Loc);
if isempty(ChannelMat) || ((~isfield(ChannelMat, 'Channel') || isempty(ChannelMat.Channel)) && ~isHeadPoints)
    disp('BST> Warning: No channel information was read from the file.');
    bst_progress('stop');
    return
end
% Are the SCS coordinates defined for this file?
isScsDefined = isfield(ChannelMat, 'SCS') && all(isfield(ChannelMat.SCS, {'NAS','LPA','RPA'})) && (length(ChannelMat.SCS.NAS) == 3) && (length(ChannelMat.SCS.LPA) == 3) && (length(ChannelMat.SCS.RPA) == 3);
% Use world coordinates by defaults for some specific file formats
if isempty(isApplyVox2ras)
    if ismember(FileFormat, {'ASCII_XYZ_WORLD', 'ASCII_NXYZ_WORLD', 'ASCII_XYZN_WORLD', 'SIMNIBS'})
        isApplyVox2ras = 1;   % Use the current vox2ras matrix in the MRI file
    elseif ismember(FileFormat, {'ASCII_XYZ', 'ASCII_NXYZ', 'ASCII_XYZN', 'BIDS-ALS-MM', 'BIDS-ALS-CM', 'BIDS-ALS-M'})
        isApplyVox2ras = 0;   % Disable vox2ras for ASCII formats that are explicitly in SCS
    elseif ismember(FileFormat, {'BIDS-SCANRAS-MM', 'BIDS-SCANRAS-CM', 'BIDS-SCANRAS-M'})
        isApplyVox2ras = 2;   % Use the vox2ras matrix AND reverts the registration done in Brainstorm, to match the original file
    end
end

%% ===== CHECK DISTANCE UNITS =====
if ~isempty(FileUnits) && ~isequal(isFixUnits, 0)
    if isempty(isFixUnits)
        isConfirmFix = 1;
    else
        isConfirmFix = 0;
    end
    ChannelMat = channel_fixunits(ChannelMat, FileUnits, isConfirmFix);
end

%% ===== MNI TRANSFORMATION =====
if ismember(FileFormat, {'ASCII_XYZ_MNI', 'ASCII_NXYZ_MNI', 'ASCII_XYZN_MNI', 'INTRANAT_MNI', 'BIDS-MNI-MM', 'BIDS-MNI-CM', 'BIDS-MNI-M'})
    % Warning for multiple studies
    if (length(iStudies) > 1)
        warning(['WARNING: When importing MNI positions for multiple subjects: the MNI transformation from the first subject is used for all of them.' 10 ...
                 'Please consider importing your subjects seprately.']);
    end
    % If we know the destination study: convert from MNI to SCS coordinates
    if ~isempty(iStudies)
        % Get the subject for the first study
        sStudy = bst_get('Study', iStudies(1));
        sSubject = bst_get('Subject', sStudy.BrainStormSubject);
        % Get the subject's MRI
        if isempty(sSubject.Anatomy) || isempty(sSubject.Anatomy(1).FileName)
            error('You need the subject anatomy in order to load sensor positions in MNI coordinates.');
        end
        % Load the MRI
        MriFile = file_fullpath(sSubject.Anatomy(1).FileName);
        sMri = in_mri_bst(MriFile);
        if ~isfield(sMri, 'SCS') || ~isfield(sMri.SCS, 'R') || isempty(sMri.SCS.R) || ~isfield(sMri, 'NCS') || ((~isfield(sMri.NCS, 'R') || isempty(sMri.NCS.R)) && (~isfield(sMri.NCS, 'y') || isempty(sMri.NCS.y)))
            error(['The SCS and MNI transformations must be defined for this subject' 10 'in order to load sensor positions in MNI coordinates.']);
        end
        % Convert all the coordinates: MNI => SCS
        fcnTransf = @(Loc)cs_convert(sMri, 'mni', 'scs', Loc')';
        AllChannelMats = channel_apply_transf(ChannelMat, fcnTransf, [], 1);
        ChannelMat = AllChannelMats{1};
    end
    % Do not convert the positions to SCS
    isAlignScs = 0;

%% ===== ACPC TRANSFORMATION =====
elseif ismember(FileFormat, {'BIDS-ACPC-MM', 'BIDS-ACPC-CM', 'BIDS-ACPC-M'})
    % Warning for multiple studies
    if (length(iStudies) > 1)
        warning(['WARNING: When importing ACPC positions for multiple subjects: the ACPC transformation from the first subject is used for all of them.' 10 ...
                 'Please consider importing your subjects seprately.']);
    end
    % If we know the destination study: convert from ACPC to SCS coordinates
    if ~isempty(iStudies)
        % Get the subject for the first study
        sStudy = bst_get('Study', iStudies(1));
        sSubject = bst_get('Subject', sStudy.BrainStormSubject);
        % Get the subject's MRI
        if isempty(sSubject.Anatomy) || isempty(sSubject.Anatomy(1).FileName)
            error('You need the subject anatomy in order to load sensor positions in ACPC coordinates.');
        end
        % Load the MRI
        MriFile = file_fullpath(sSubject.Anatomy(1).FileName);
        sMri = in_mri_bst(MriFile);
        if ~isfield(sMri, 'SCS') || ~isfield(sMri.SCS, 'R') || isempty(sMri.SCS.R) || ~isfield(sMri.NCS, 'AC') || isempty(sMri.NCS.AC) || ~isfield(sMri.NCS, 'PC') || isempty(sMri.NCS.PC) || ~isfield(sMri.NCS, 'IH') || isempty(sMri.NCS.IH)
            error(['All fiducials must be defined for this subject (NAS,LPA,RPA,AC,PC,IH)' 10 'in order to load sensor positions in ACPC coordinates.']);
        end
        % Convert all the coordinates: ACPC => SCS
        fcnTransf = @(Loc)cs_convert(sMri, 'acpc', 'scs', Loc')';
        AllChannelMats = channel_apply_transf(ChannelMat, fcnTransf, [], 1);
        ChannelMat = AllChannelMats{1};
    end
    % Do not convert the positions to SCS
    isAlignScs = 0;

%% ===== CAPTRAK TRANSFORMATION =====
elseif ismember(FileFormat, {'BIDS-CAPTRAK-MM', 'BIDS-CAPTRAK-CM', 'BIDS-CAPTRAK-M'})
    % Warning for multiple studies
    if (length(iStudies) > 1)
        warning(['WARNING: When importing CapTrak positions for multiple subjects: the CapTrak transformation from the first subject is used for all of them.' 10 ...
                 'Please consider importing your subjects seprately.']);
    end
    % If we know the destination study: convert from CapTrak to SCS coordinates
    if ~isempty(iStudies)
        % Get the subject for the first study
        sStudy = bst_get('Study', iStudies(1));
        sSubject = bst_get('Subject', sStudy.BrainStormSubject);
        % Get the subject's MRI
        if isempty(sSubject.Anatomy) || isempty(sSubject.Anatomy(1).FileName)
            error('You need the subject anatomy in order to load sensor positions in CapTrak coordinates.');
        end
        % Load the MRI
        MriFile = file_fullpath(sSubject.Anatomy(1).FileName);
        sMri = in_mri_bst(MriFile);
        if ~isfield(sMri, 'SCS') || ~isfield(sMri.SCS, 'R') || isempty(sMri.SCS.R) || ~isfield(sMri.SCS, 'NAS') || isempty(sMri.SCS.NAS) || ~isfield(sMri.SCS, 'LPA') || isempty(sMri.SCS.LPA) || ~isfield(sMri.SCS, 'RPA') || isempty(sMri.SCS.RPA)
            error(['All fiducials must be defined for this subject (NAS,LPA,RPA)' 10 'in order to load sensor positions in CapTrak coordinates.']);
        end
        % Convert all the coordinates: CapTrak => SCS
        fcnTransf = @(Loc)cs_convert(sMri, 'captrak', 'scs', Loc')';
        AllChannelMats = channel_apply_transf(ChannelMat, fcnTransf, [], 1);
        ChannelMat = AllChannelMats{1};
    end
    % Do not convert the positions to SCS
    isAlignScs = 0;

%% ===== MRI/NII TRANSFORMATION =====
% If the SCS coordinates are not defined (NAS/LPA/RPA fiducials), try to use the MRI=>subject transformation available in the MRI (eg. NIfTI sform/qform)
% Only available if there is one study in output
elseif ~isScsDefined && ~isequal(isApplyVox2ras, 0) && ~isempty(iStudies)
    % Get the folders
    sStudies = bst_get('Study', iStudies);
    if (length(sStudies) > 1) && ~all(strcmpi(sStudies(1).BrainStormSubject, {sStudies.BrainStormSubject}))
        warning(['WARNING: When importing sensor positions for multiple subjects: the SCS transformation from the first subject is used for all of them.' 10 ...
                 'Please consider importing your subjects seprately.']);
    end
    % Get subject for first file only
    sSubject = bst_get('Subject', sStudies(1).BrainStormSubject);

    % SIMNIBS: Disable auto-alignment based on fiducials
    if strcmpi(FileFormat, 'SIMNIBS')
        isAlignScs = 0;
    else
        isAlignScs = 1;
    end
    % If there is a MRI for this subject
    if ~isempty(sSubject.Anatomy) && ~isempty(sSubject.Anatomy(1).FileName)
        % Get the reference MRI (specified in input, or selected in the database)
        if isempty(RefMriFile) || ~file_exist(file_fullpath(RefMriFile))
            % If there are multiple MRIs
            if (length(sSubject.Anatomy) > 1)
                % Consider only the volumes that are not volume atlases
                iNoAtlas = find(cellfun(@(c)isempty(strfind(c, '_volatlas')), {sSubject.Anatomy.FileName}));
                if (length(iNoAtlas) == 1)
                    iMri = iNoAtlas;
                % Interactive: Ask which MRI volume to use
                elseif isInteractive
                    mriComment = java_dialog('combo', '<HTML>Select the reference MRI:<BR><BR>', 'Import as MRI scanner coordinates', [], {sSubject.Anatomy(iNoAtlas).Comment});
                    if isempty(mriComment)
                        bst_progress('stop');
                        return
                    end
                    iMri = iNoAtlas(find(strcmp({sSubject.Anatomy(iNoAtlas).Comment}, mriComment), 1));
                % Non-interactive: Use the default MRI
                else
                    iMri = sSubject.iAnatomy;
                end
            else
                iMri = 1;
            end
            RefMriFile = file_fullpath(sSubject.Anatomy(iMri).FileName);
        end
        % Load the reference MRI (which contains the vox2mri transformation that should be used to interpret the coordinates)
        sMri = load(file_fullpath(RefMriFile), 'InitTransf', 'SCS', 'Voxsize');
        % If there is a valid transformation
        if isfield(sMri, 'InitTransf') && ~isempty(sMri.InitTransf) && ismember('vox2ras', sMri.InitTransf(:,1))
            % Ask user if necessary
            if isempty(isApplyVox2ras)
                isApplyVox2ras = java_dialog('confirm', [...
                    'There is a transformation to subject coordinates available in the MRI.' 10 ...
                    'Would you like to use it to align the sensors with the MRI?' 10 10 ...
                    'Answer NO to use the NAS/LPA/RPA fiducials from the input file.' 10], 'Apply MRI transformation');
            end
            % Apply transformation
            if isApplyVox2ras
                % Get the transformation WORLD=>MRI (in meters)
                Transf = cs_convert(sMri, 'world', 'mri');
                % Applies the coregistration transformation to the sensor positions, if requested
                if (isApplyVox2ras == 2)
                    iTransfReg = find(strcmpi(sMri.InitTransf(:,1), 'reg'), 1);
                    if ~isempty(iTransfReg)
                        Transf = Transf * sMri.InitTransf{iTransfReg,2};
                    end
                end
                % Add the transformation MRI=>SCS
                if isfield(sMri,'SCS') && isfield(sMri.SCS,'R') && ~isempty(sMri.SCS.R) && isfield(sMri.SCS,'T') && ~isempty(sMri.SCS.T)
                    Transf = [sMri.SCS.R, sMri.SCS.T./1000; 0 0 0 1] * Transf;
                else
                    error(['The SCS coordinates are not defined for this subject, the sensors will not be aligned on the anatomy. ' 10 'Consider defining the NAS/LPA/RPA fiducials before importing the sensors positions.']);
                end
                % Convert all the coordinates
                AllChannelMats = channel_apply_transf(ChannelMat, Transf, [], 1);
                ChannelMat = AllChannelMats{1};
                % Disable alignment based on SCS fiducials
                isAlignScs = 0;
            end
        end
    end

elseif isfield(ChannelMat, 'TransfMegLabels') && iscell(ChannelMat.TransfMegLabels) && ismember('Native=>Brainstorm/CTF', ChannelMat.TransfMegLabels)
    % No need to duplicate this transformation if it was previously
    % computed, e.g. in in_channel_ctf. (It would be identity the second
    % time.)
    isAlignScs = 0;
else
    isAlignScs = 1;
end

%% ===== DETECT CHANNEL TYPES =====
% Remove fiducials (expect for BIDS files)
isRemoveFid = isempty(strfind(FileFormat, 'BIDS-'));
% Detect auxiliary EEG channels + align channel
ChannelMat = channel_detect_type(ChannelMat, isAlignScs, isRemoveFid);

%% ===== APPLY NEW CHANNEL FILE =====
% If some studies were defined
if isSave && ~isempty(iStudies)
    if isempty(ChannelAlign)
        iMEG = good_channel(ChannelMat.Channel, [], 'MEG');
        ChannelAlign = ~isempty(iMEG);
    end
    % History: Import channel file
    ChannelMat = bst_history('add', ChannelMat, 'import', ['Import from: ' ChannelFile ' (Format: ' FileFormat ')']);
    % Add channel file to all the target studies
    for i = 1:length(iStudies)
        ChannelFile = db_set_channel(iStudies(i), ChannelMat, ChannelReplace, ChannelAlign);
    end
    % Returned value
    Output = ChannelFile;
else
    Output = ChannelMat;
end

% Progress bar
if ~isProgressBar
    bst_progress('stop');
end

author: brainstorm-tools


text:
function msm_to_mm_test17 ( )

%*****************************************************************************80
%
%% MSM_TO_MM_TEST17 tests MSM_TO_MM_COORDINATE_INTEGER_SYMMETRIC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MSM_TO_MM_TEST17\n' );
  fprintf ( 1, '  Convert an MSM to MM coordinate integer symmetric format.\n' );

  output_filename = 'msm_to_mm_test17.mm';

  a = i4mat_indicator ( 4, 4 );
  a = a + a';
%
%  Have MSM_TO_MM write the matrix to a file.
%
  msm_to_mm ( output_filename, a, 'coordinate', 'integer', 'symmetric' );

  return
end


text:
clear all; close all; clc;
%% Semi-Supervised Generative Adversarial Network
%% Load Data
load('mnistAll.mat')
trainX = preprocess(mnist.train_images);
trainY = mnist.train_labels;
testX = preprocess(mnist.test_images);
testY = mnist.test_labels;
%% Settings
settings.latentDim = 100; settings.num_classes = 10;
settings.batch_size = 32; settings.image_size = [28,28,1];
settings.lr = 0.0002; settings.beta1 = 0.5;
settings.loss_weights = [.5,.5];
settings.beta2 = 0.999; settings.maxepochs = 50;

%% Initialization
%% Generator
paramsGen.FCW1 = dlarray(initializeGaussian([128*7*7,...
    settings.latentDim]));
paramsGen.FCb1 = dlarray(zeros(128*7*7,1,'single'));
paramsGen.BNo1 = dlarray(zeros(128,1,'single'));
paramsGen.BNs1 = dlarray(ones(128,1,'single'));
paramsGen.TCW1 = dlarray(initializeGaussian([3,3,128,128]));
paramsGen.TCb1 = dlarray(zeros(128,1,'single'));
paramsGen.BNo2 = dlarray(zeros(128,1,'single'));
paramsGen.BNs2 = dlarray(ones(128,1,'single'));
paramsGen.TCW2 = dlarray(initializeGaussian([3,3,64,128]));
paramsGen.TCb2 = dlarray(zeros(64,1,'single'));
paramsGen.BNo3 = dlarray(zeros(64,1,'single'));
paramsGen.BNs3 = dlarray(ones(64,1,'single'));
paramsGen.CNW1 = dlarray(initializeGaussian([3,3,64,1]));
paramsGen.CNb1 = dlarray(zeros(1,1,'single'));
stGen.BN1 = []; stGen.BN2 = []; stGen.BN3 = [];

%% Discriminator
paramsDis.CNW1 = dlarray(initializeGaussian([3,3,1,32]));
paramsDis.CNb1 = dlarray(zeros(32,1,'single'));
paramsDis.CNW2 = dlarray(initializeGaussian([3,3,32,64]));
paramsDis.CNb2 = dlarray(zeros(64,1,'single'));
paramsDis.BNo1 = dlarray(zeros(64,1,'single'));
paramsDis.BNs1 = dlarray(ones(64,1,'single'));
paramsDis.CNW3 = dlarray(initializeGaussian([3,3,64,128]));
paramsDis.CNb3 = dlarray(zeros(128,1,'single'));
paramsDis.BNo2 = dlarray(zeros(128,1,'single'));
paramsDis.BNs2 = dlarray(ones(128,1,'single'));
paramsDis.CNW4 = dlarray(initializeGaussian([3,3,128,256]));
paramsDis.CNb4 = dlarray(zeros(256,1,'single'));
paramsDis.FCW1 = dlarray(initializeGaussian([settings.num_classes+2,256*4*4]));
paramsDis.FCb1 = dlarray(zeros(settings.num_classes+2,1,'single'));
stDis.BN1 = []; stDis.BN2 = []; stDis.BN3 = [];

% average Gradient and average Gradient squared holders
avgG.Dis = []; avgGS.Dis = []; avgG.Gen = []; avgGS.Gen = [];
%% Train
numIterations = floor(size(trainX,4)/settings.batch_size);
out = false; epoch = 0; global_iter = 0;
while ~out
    tic;
    shuffleid = randperm(size(trainX,4));
    trainXshuffle = trainX(:,:,:,shuffleid);
    trainYshuffle = trainY(shuffleid);
    fprintf('Epoch %d\n',epoch)
    for i=1:numIterations
        global_iter = global_iter+1;
        noise = gpdl(randn([settings.latentDim,...
            settings.batch_size]),'CB');
        idx = (i-1)*settings.batch_size+1:i*settings.batch_size;
        XBatch=gpdl(single(trainXshuffle(:,:,:,idx)),'SSCB');
        YBatch=gpdl(single(trainYshuffle(idx)),'B');

        [GradGen,GradDis,stGen,stDis] = ...
                dlfeval(@modelGradients,XBatch,YBatch,noise,...
                paramsGen,paramsDis,stGen,stDis);

        % Update Discriminator network parameters
        [paramsDis,avgG.Dis,avgGS.Dis] = ...
            adamupdate(paramsDis, GradDis, ...
            avgG.Dis, avgGS.Dis, global_iter, ...
            settings.lr, settings.beta1, settings.beta2);

        % Update Generator network parameters
        [paramsGen,avgG.Gen,avgGS.Gen] = ...
            adamupdate(paramsGen, GradGen, ...
            avgG.Gen, avgGS.Gen, global_iter, ...
            settings.lr, settings.beta1, settings.beta2);

        if i==1 || rem(i,20)==0
            progressplot(paramsGen,stGen,settings);
        end

    end

    elapsedTime = toc;
    disp("Epoch "+epoch+". Time taken for epoch = "+elapsedTime + "s")
    epoch = epoch+1;
    if epoch == settings.maxepochs
        out = true;
    end
end
%% Helper Functions
%% one hot encoding
function ohe = onehotencoding(labels,numLabels)
numBatch = length(labels);
ohe = zeros(numLabels,numBatch);

for i = 1:numBatch
    ohe(labels(i)+1,i)=1;
end
ohe(1:10,:) = 5/8*ohe(1:10,:);
ohe(end,:) = 1/16*ohe(end,:);

end
%% preprocess
function x = preprocess(x)
x = double(x)/255;
x = (x-.5)/.5;
x = reshape(x,28,28,1,[]);
end
%% extract data
function x = gatext(x)
x = gather(extractdata(x));
end
%% gpu dl array wrapper
function dlx = gpdl(x,labels)
dlx = gpuArray(dlarray(x,labels));
end
%% Weight initialization
function parameter = initializeGaussian(parameterSize,sigma)
if nargin < 2
    sigma = 0.05;
end
parameter = randn(parameterSize, 'single') .* sigma;
end
%% Generator
function [dly,st] = Generator(dlx,params,st)
% fully connected
dly = fullyconnect(dlx,params.FCW1,params.FCb1);
dly = relu(dly);
% transposed convolution
dly = gpdl(reshape(dly,7,7,128,[]),'SSCB');
if isempty(st.BN1)
    [dly,st.BN1.mu,st.BN1.sig] = batchnorm(dly,...
        params.BNo1,params.BNs1,'MeanDecay',0.8);
else
    [dly,st.BN1.mu,st.BN1.sig] = batchnorm(dly,params.BNo1,...
        params.BNs1,st.BN1.mu,st.BN1.sig,...
        'MeanDecay',.8);
end
dly = dltranspconv(dly,params.TCW1,params.TCb1,...
    'Stride',2,'Cropping','same');
dly = relu(dly);
if isempty(st.BN2)
    [dly,st.BN2.mu,st.BN2.sig] = batchnorm(dly,...
        params.BNo2,params.BNs2,'MeanDecay',0.8);
else
    [dly,st.BN2.mu,st.BN2.sig] = batchnorm(dly,params.BNo2,...
        params.BNs2,st.BN2.mu,st.BN2.sig,...
        'MeanDecay',.8);
end

dly = dltranspconv(dly,params.TCW2,params.TCb2,...
    'Stride',2,'Cropping','same');
dly = relu(dly);
if isempty(st.BN3)
    [dly,st.BN3.mu,st.BN3.sig] = batchnorm(dly,...
        params.BNo3,params.BNs3,'MeanDecay',0.8);
else
    [dly,st.BN3.mu,st.BN3.sig] = batchnorm(dly,params.BNo3,...
        params.BNs3,st.BN3.mu,st.BN3.sig,...
        'MeanDecay',.8);
end

dly = dlconv(dly,params.CNW1,params.CNb1,...
            'Padding','same');
% tanh
dly = tanh(dly);
end
%% Discriminator
function [dly,st] = Discriminator(dlx,params,st)
% convolution
%1
dly = dlconv(dlx,params.CNW1,params.CNb1,...
            'Stride',2,'Padding','same');
dly = leakyrelu(dly,0.2);
dly = dropout(dly,.25);
%2
dly = dlconv(dly,params.CNW2,params.CNb2,...
            'Stride',2,'Padding','same');
dly = leakyrelu(dly,0.2);
dly = dropout(dly,.25);
if isempty(st.BN1)
    [dly,st.BN1.mu,st.BN1.sig] = batchnorm(dly,...
        params.BNo1,params.BNs1,'MeanDecay',0.8);
else
    [dly,st.BN1.mu,st.BN1.sig] = batchnorm(dly,params.BNo1,...
        params.BNs1,st.BN1.mu,st.BN1.sig,...
        'MeanDecay',0.8);
end
%3
dly = dlconv(dly,params.CNW3,params.CNb3,...
            'Stride',2,'Padding','same');
dly = leakyrelu(dly,0.2);
dly = dropout(dly,.25);
if isempty(st.BN2)
    [dly,st.BN2.mu,st.BN2.sig] = batchnorm(dly,...
        params.BNo2,params.BNs2,'MeanDecay',0.8);
else
    [dly,st.BN2.mu,st.BN2.sig] = batchnorm(dly,params.BNo2,...
        params.BNs2,st.BN2.mu,st.BN2.sig,...
        'MeanDecay',0.8);
end
%4
dly = dlconv(dly,params.CNW4,params.CNb4,...
            'Stride',1,'Padding','same');
dly = leakyrelu(dly,0.2);
dly = dropout(dly,.25);

% Fully connected
dly = gpdl(reshape(dly,4*4*256,[]),'CB');
dly = fullyconnect(dly,params.FCW1,params.FCb1);
% sigmoid
dly(1,:) = sigmoid(dly(1,:));
% softmax
dly(2:end,:) = softmax(dly(2:end,:));
end
%% modelGradients
function [GradGen,GradDis,stGen,stDis]=modelGradients(x,y,z,paramsGen,...
    paramsDis,stGen,stDis)
[fake_images,stGen] = Generator(z,paramsGen,stGen);
d_output_real = Discriminator(x,paramsDis,stDis);
[d_output_fake,stDis] = Discriminator(fake_images,paramsDis,stDis);
labels_real = onehotencoding(y,11);
labels_fake = onehotencoding(10*ones(length(y),1),11);

% Loss due to true or not
dlossreal=-mean(log(d_output_real(1,:)+eps)+...
    labels_real.*log(d_output_real(2:end,:)+eps),'all');
dlossfake=-mean(log(1-d_output_fake(1,:)+eps)+...
    labels_fake.*log(d_output_fake(2:end,:)+eps),'all');
d_loss = .5*(dlossreal+dlossfake);
% g_loss=-mean(log(d_output_fake(1,:)+eps)+...
%     labels_real.*log(d_output_fake(2:end,:)+eps),'all');
g_loss=-mean(log(d_output_fake(1,:)+eps),'all');

% For each network, calculate the gradients with respect to the loss.
GradGen = dlgradient(g_loss,paramsGen,'RetainData',true);
GradDis = dlgradient(d_loss,paramsDis);
end
%% progressplot
function progressplot(paramsGen,stGen,settings)
r = 5; c = 5;
noise = gpdl(randn([settings.latentDim,r*c]),'CB');
gen_imgs = Generator(noise,paramsGen,stGen);
gen_imgs = reshape(gen_imgs,28,28,[]);

fig = gcf;
if ~isempty(fig.Children)
    delete(fig.Children)
end

I = imtile(gatext(gen_imgs));
I = rescale(I);
imagesc(I)
title("Generated Images")
colormap gray

drawnow;
end
%% dropout
function dly = dropout(dlx,p)
if nargin < 2
    p = .3;
end
[n,d] = rat(p);
mask = randi([1,d],size(dlx));
mask(mask<=n)=0;
mask(mask>n)=1;
dly = dlx.*mask;

end

author: zcemycl

repo: Matlab-GAN

save_path: github-repos/MATLAB/zcemycl-Matlab-GAN


text:
%%
close all; clc;
addpath('./scripts');
col = {[0 114 178],[0 158 115], [213 94 0],[230 159 0],...
    [86 180 233], [204 121 167], [64 224 208], [240 228 66]};
figure('papersize', [14, 7]);
init_fig;

author: zhoupc


text:
%% housekeeping
clear
close all
clc()

%% Create dataset
clc

do_plot=true;

scale=100;

[db,varlist0]=create_dataset(scale,do_plot);

varlist=fieldnames(varlist0);
%% Choose a model type: see cell "create the structural VAR model" below
model_type=4;

%% set up the restrictions
close()

% create restrictions on parameters as well as markov chains
%------------------------------------------------------------
switch model_type
    case 0
        % constant-parameter model
        [lin_restr,nonlin_restr,markov_chains]=create_restrictions_and_markov_chains0();
    case 1
        % Coefficients are switching regimes across all equations
        % (synchronized case)
        [lin_restr,nonlin_restr,markov_chains]=create_restrictions_and_markov_chains1();
    case 2
        % Coefficients and variances have different chains, different
        % regimes, and different durations
        [lin_restr,nonlin_restr,markov_chains]=create_restrictions_and_markov_chains2();
    case 3
        % Only coefficients in monetary policy equation are changing
        [lin_restr,nonlin_restr,markov_chains]=create_restrictions_and_markov_chains3();
    case 4
        % Only variance in monetary policy equation is changing
        [lin_restr,nonlin_restr,markov_chains]=create_restrictions_and_markov_chains4();
    case 5
        % Both coefficients and variances in monetary policy equation
        % change with two independent Markov processes
        [lin_restr,nonlin_restr,markov_chains]=create_restrictions_and_markov_chains5();
    case 6 % ok
        % Only variances in ALL three equations switch
        [lin_restr,nonlin_restr,markov_chains]=create_restrictions_and_markov_chains6();
    otherwise
        error('the coded model types are 0, 1, 2, 3, 4 and 6')
end

%% Create the VAR

nlags=2;

exog={};

panel=[];

constant=true;

% first we create a template structure
% ------------------------------------
sv0=svar(varlist,exog,nlags,constant,panel,markov_chains);

%% set priors

prior=svar.prior_template();

prior.type='sz';

% prior.L1=0.1/2;
%
% prior.coefprior=0.5;

is_prior=~true;

%% Find posterior mode
clc

sv=sv0;

if is_prior

    sv=set(sv,'prior',prior);

end

sv=estimate(sv,{'1960Q1','2015Q2'},'data',db,...
    'linear_restrictions',[lin_restr;nonlin_restr]);

%% Printing estimates

print_structural_form(sv)

%% Printing solution
clc

print_solution(sv)

%% plot smoothed state and regime probabilities

close all

plot_probabilities(sv)

%% plots probabilities against data
close all

plot_data_against_probabilities(sv,'state')

%% Impulse responses

myirfs=irf(sv);

%% Posterior sampling

%% Marginal data density

%% Out-of sample forecasts

%% Conditional forecast

author: jmaih


text:
function CPD = update_CPT(CPD)
% Compute the big CPT for an HHMM Q node (including F parents) given internal transprob and startprob
% function CPD = update_CPT(CPD)

Qsz = CPD.Qsz;
Qpsz = CPD.Qpsz;

if ~isempty(CPD.Fbelow_ndx)
  if ~isempty(CPD.Fself_ndx) % general case
    % Fb(t-1) Fself(t-1)  P(Q(t)=j| Q(t-1)=i, Qps(t)=k)
    % ------------------------------------------------------
    % 1        1         delta(i,j)
    % 2        1         transprob(i,k,j)
    % 1        2         impossible
    % 2        2         startprob(k,j)
    CPT = zeros(Qsz, 2, 2, Qpsz, Qsz);
    I = repmat(eye(Qsz), [1 1 Qpsz]); % i,j,k
    I = permute(I, [1 3 2]); % i,k,j
    CPT(:, 1, 1, :, :) = I;
    CPT(:, 2, 1, :, :) = CPD.transprob;
    CPT(:, 1, 2, :, :) = I;
    CPT(:, 2, 2, :, :) = repmat(reshape(CPD.startprob, [1 Qpsz Qsz]), [Qsz 1 1]); % replicate over i
  else % no F from self, hence no startprob
    % Fb(t-1) P(Q(t)=j| Q(t-1)=i, Qps(t)=k)
    % ------------------------------------------------------
    % 1       delta(i,j)
    % 2       transprob(i,k,j)

    nps = length(CPD.dom_sz)-1; % num parents
    CPT = 0*myones(CPD.dom_sz);
    %CPT = zeros(Qsz, 2, Qpsz, Qsz); % assumes CPT(Q(t-1), F(t-1), Qps, Q(t))
    % but a member of Qps may preceed Q(t-1) or F(t-1) in the ordering

    I = repmat(eye(Qsz), [1 1 Qpsz]); % i,j,k
    I = permute(I, [1 3 2]); % i,k,j

    % the following fails if there is a member of Qps with a lower
    % number than F
    %CPT(:, 1, :, :) = I;
    %CPT(:, 2, :, :) = CPD.transprob;

    ndx = mk_multi_index(nps+1, CPD.Fbelow_ndx, 1);
    CPT(ndx{:}) = I;
    ndx = mk_multi_index(nps+1, CPD.Fbelow_ndx, 2);
    CPT(ndx{:}) = CPD.transprob;
    keyboard
  end
else % no F signal from below
  if ~isempty(CPD.Fself_ndx)
    % Q(t-1), Fself(t-1), Qps, Q(t)

    % if condition start on previous concrete state (as in map learning),
    % CPT(:, 1, :, :, :) = CPD.transprob(Q(t-1), Qps, Q(t))
    % CPT(:, 2, :, :, :) = CPD.startprob(Q(t-1), Qps, Q(t))

    % Fself(t-1)  P(Q(t-1)=i, Qps(t)=k -> Q(t)=j)
    % ------------------------------------------------------
    % 1         transprob(i,k,j)
    % 2         startprob(k,j)
    CPT = zeros(Qsz, 2, Qpsz, Qsz);
    I = repmat(eye(Qsz), [1 1 Qpsz]); % i,j,k
    I = permute(I, [1 3 2]); % i,k,j
    CPT(:, 1, :, :) = CPD.transprob;
    if CPD.fullstartprob
      CPT(:, 2, :, :) = CPD.startprob;
    else
      CPT(:, 2, :, :) = repmat(reshape(CPD.startprob, [1 Qpsz Qsz]), [Qsz 1 1]); % replicate over i
    end
    else % no F from self
    error('An hhmmQ node without any F parents is just a tabular_CPD')
  end
end

CPD = set_fields(CPD, 'CPT', CPT);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function f1 = perform_curve_subdivision(f, nsub, options)

% perform_curve_subdivision - perform subdivision
%
%   f1 = perform_curve_subdivision(f, nsub, options);
%
%   Copyright (c) 2007 Gabriel Peyre

options.null = 0;
if nsub<=0
    f1 = f; return;
end
if size(f,1)>size(f,2)
    f = f';
end
if size(f,1)>1
    for k=1:size(f,1)
        f1(k,:) = perform_curve_subdivision(f(k,:), nsub, options);
    end
    return;
end

h = getoptions(options, 'h', [1 4 6 4 1]);
h = 2*h/sum(h(:));

if not(isfield(options, 'bound'))
    options.bound = 'per';
end

f = perform_curve_subdivision(f, nsub-1, options);

f = f(:)';
n = length(f);
f1 = zeros(1,n*2);
f1(1:2:end) = f;
f1 = perform_convolution(f1, h, options);
f1 = f1(:)';

author: gpeyre

repo: matlab-toolboxes


text:
%
%   FILE NAME:    HMR2Curation.m
%
%   PURPOSE: Curating HMR2 database model toward HMR3
%

% Load the original Matlab file of HMR2 database
load('HMRdatabase2_00.mat');

% Import HMR2 reactions from the Excel file 'RXNS' sheet that leave out
% the following five rows that represent reaction classes
% 7709: 'Exchange reactions';
% 8171: 'Fake reactions';
% 8176: 'Biomass reactions';
% 8177: 'HMR_biomass_Renalcancer';
% 8183: 'Included for connectivity for INIT'
T=readtable('HMRdatabase2_00.xlsx','Sheet','RXNS','ReadVariableNames',1);
HMR2=table2struct(T,'ToScalar',true);

% Adding reaction identifiers from external databases
% additional empty spaces were also removed druing this process
if isequal(ihuman.rxns,HMR2.RXNID)
	ihuman.rxnKEGGID=HMR2.KEGGID;              %KEGG
	ihuman.rxnEHMNID=HMR2.EHMNID;              %EHMN
	ihuman.rxnBiGGID=HMR2.BIGGDATABASEID;      %Recon
	ihuman.rxnHepatoNET1ID=HMR2.HEPATONET1ID;  %Hepatonet1
	ihuman.rxnREACTOMEID=HMR2.REACTOMEID;      %Reactome
	ihuman.rxnReferences=HMR2.REFERENCES;      %References
end

% Some errors were spotted and fixed:
% The reaction HMR_2190 was assocated to four REACTOM reactions.
% The first one (REACT_22270) is correect, the other three
% (REACT_22097; REACT_22133; REACT_22219) should be wrong because
% they are empty ids by searching reactom.org, and thus removed
index=find(strcmp('HMR_2190',ihuman.rxns));
ihuman.rxnREACTOMEID{index}='REACT_22270';

% The KEGG id of HMR_7709 was associated to 'R0302' that should be typo
% It was manually checked and corrected to R03027
index=find(strcmp('HMR_7709',ihuman.rxns));
ihuman.rxnKEGGID{index}='R03027';

% Save as version 2.0.1
save('HMRdatabase2_01.mat','ihuman'); %===2018-01-16

author: SysBioChalmers

repo: Human-GEM


text:
function test_ft_plot_vector

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_plot_vector

mask = zeros(1, 100)
%boundary conditions
mask(1:10) = 1;
mask(90:100) = 1;
mask(40:60) = 1;

figure
subplot(2,2,1)
ft_plot_vector(randn(1,100), 'width', 1, 'height', 1, 'hpos', 0, 'vpos', 0);
subplot(2,2,2)
ft_plot_vector(randn(1,100), 'width', 1, 'height', 1, 'hpos', 0, 'vpos', 0, 'highlight', mask);
subplot(2,2,3)
ft_plot_vector(randn(1,100), 'width', 1, 'height', 1, 'hpos', 0, 'vpos', 0, 'highlight', mask, 'highlightstyle', 'saturation');
subplot(2,2,4)
ft_plot_vector(randn(1,100), 'width', 1, 'height', 1, 'hpos', 0, 'vpos', 0, 'highlight', mask, 'highlightstyle', 'thickness');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%u2F    estimates fundamental matrix using ortogonal LS regression
%	F = u2F(u) estimates F from u using NORMU
%  F = u2F(u,'nonorm') disables normalization
%	see also NORMU, U2FA

function F = u2F (u, str)

if (nargin > 1) & strcmp(str, 'nonorm')
   donorm = 0;
else
   donorm = 1;
end

ptNum = size(u,2);

if donorm
   A1    = normu(u(1:3,:));
   A2    = normu(u(4:6,:));

   u1   = A1*u(1:3,:);
   u2   = A2*u(4:6,:);
end

for i = 1:ptNum
   Z(i,:)   = reshape(u1(:,i)*u2(:,i)',1,9);
end

M       = Z'*Z;
V       = seig(M);
F = reshape(V(:,1),3,3);

[uu,us,uv] = svd(F);
[y,i]      = min (abs(diag(us)));
us(i,i)    = 0;
F          = uu*us*uv';

if donorm
   F = A1'*F*A2;
end

F = F /norm(F,2);

author: strawlab

repo: MultiCamSelfCal


text:
%{
    Copyright (c) 2015, Philipp Krähenbühl
    All rights reserved.

    Redistribution and use in source and binary forms, with or without
    modification, are permitted provided that the following conditions are met:
        * Redistributions of source code must retain the above copyright
        notice, this list of conditions and the following disclaimer.
        * Redistributions in binary form must reproduce the above copyright
        notice, this list of conditions and the following disclaimer in the
        documentation and/or other materials provided with the distribution.
        * Neither the name of the Stanford University nor the
        names of its contributors may be used to endorse or promote products
        derived from this software without specific prior written permission.

    THIS SOFTWARE IS PROVIDED BY Philipp Krähenbühl ''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 Philipp Krähenbühl 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.
%}
addpath( genpath('../build') )

author: Cloud-CV

repo: object-proposals


text:
function fx2 = p18_fx2 ( x )

%*****************************************************************************80
%
%% P18_FX2 evaluates the second derivative of the function for problem 18.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(*), the point at which F is to be evaluated.
%
%    Output, real FX2(*), the value of the second derivative at X.
%
  fx2 = 10.0^14 * ( ...
         42 * x.^5 ...
      - 210 * x.^4 ...
      + 420 * x.^3 ...
      - 420 * x.^2 ...
      + 210 * x    ...
      -  42 );

  return
end


text:
function Mh = hmxChol(Mh)
%+========================================================================+
%|                                                                        |
%|         OPENHMX - LIBRARY FOR H-MATRIX COMPRESSION AND ALGEBRA         |
%|           openHmx is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : hmxChol.m                                     |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
%| ( === ) |   SYNOPSIS   : Cholesky factorization of H-Matrix            |
%|  `---'  |                                                              |
%+========================================================================+

% Check invertibility
if (size(Mh,1) ~= size(Mh,2)) || ~isequal(Mh.pos{1},Mh.pos{2})
    error('hmxChol.m : matrix must be invertible.')
end

%%% H-Matrix (recursion)
if (Mh.typ == 0)
    % U11 -> M11
    Mh.chd{1} = hmxChol(Mh.chd{1});

    % U12 -> U11' \ M12
    Mh.chd{2} = Mh.chd{1}'\Mh.chd{2};

    % U21 -> 0
    Mh.chd{3} = zeros(Mh.chd{3});

    % M22 -> M22 - U12'*U12
    Mh.chd{4} = Mh.chd{4} - Mh.chd{2}' * Mh.chd{2};
%     Mh.chd{4} = plusmtimes(Mh.chd{4},-1,Mh.chd{2}',Mh.chd{2});

    % U22 -> M22
    Mh.chd{4} = hmxChol(Mh.chd{4});

    % Fusion
    Mh = hmxFusion(Mh);

%%% Compressed leaf
elseif (Mh.typ == 1)
    error('hmxChol : unavailable case')

%%% Full leaf
elseif (Mh.typ == 2)
    Mh.dat = chol(Mh.dat);

%%% Unknown type
else
    error('hmxChol.m : unavailable case')
end
end


text:
% Map spike timestamps to some larger time vector
%
% This function return indices of ts that correspond to spikes values.
% The idea is that ts a long signal and spikes are found within it.
% One application is for example to match spikes to LFP signal.
%
%  USAGE
%   spkInd = general.spikes2time(spikes, ts)
%   spikes      Vector of spike timestamps. Units of spikes should be the
%               same as units of ts.
%   ts          Time signal. Vector or matrix. If matrix is provided,
%               the first column is taken as ts.
%   spkInd      Vector of indices in ts that correspond to spike values.
%
function spkInd = spikes2time(spikes, ts)
    if size(ts, 1) > 1 && size(ts, 2) > 1
        ts = ts(:, 1);
    else
        ts = ts(:);
    end
    % make sure it is a column
    spikes = spikes(:);

    spkInd = knnsearch(ts, spikes);
end

author: kavli-ntnu


text:
function inc = crswap ( a, c, c_size, m, k, n, critvl, i1, i2, c1, c2, ...
  iswitch )

%*****************************************************************************80
%
%% CRSWAP determines the effect of swapping two objects.
%
%  Discussion:
%
%    This computation is very inefficient.  It is only set up so that we
%    can compare algorithm ASA 113 to the K-means algorithms ASA 058 and
%    ASA 136.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Colin Banfield, LC Bassill,
%    Algorithm AS 113:
%    A transfer for non-hierarchichal classification,
%    Applied Statistics,
%    Volume 26, Number 2, 1977, pages 206-210.
%
%  Parameters:
%
%    Input, real  A(M,N), the data values.  There are M objects,
%    each having spatial dimension N.
%
%    Input, integer C(M), the classification of each object.
%
%    Input, integer  C_SIZE(K), the number of objects in each class.
%
%    Input, integer M, the number of objects.
%
%    Input, integer K, the number of classes.
%
%    Input, integer N, the number of spatial dimensions, or variates,
%    of the objects.
%
%    Input, real CRITVL, the current value of the criterion.
%
%    Input, integer I1, I2, the objects to be swapped.
%
%    Input, integer C1, C2, the current classes of objects I1 and I2.
%
%    Input, integer ISWITCH:
%    1, indicates that I1 and I2 should be temporarily swapped, the
%       change in CRITVL should be computed, and then I1 and I2 restored.
%    2, indicates that I1 and I2 will be swapped.
%
%    Output, real INC, the change to CRITVL that would occur if I1 and
%    I2 were swapped.  This is only computed for ISWITCH = 1.
%
  if ( iswitch == 2 )
    inc = 0;
    return
  end
%
%  Move object I1 from class C1 to class C2.
%  Move object I2 from class C2 to class C1.
%
  c(i1) = c2;
  c(i2) = c1;
%
%  Define the critical value as the sum of the squares of the distances
%  of the points to their cluster center.
%
  for i = 1 : k
    c_size(i) = 0;
    for j = 1 : n
      c_center(i,j) = 0.0;
    end
  end

  for i = 1 : m
    ci = c(i);
    c_size(ci) = c_size(ci) + 1;
    for j = 1 : n
      c_center(ci,j) = c_center(ci,j) + a(i,j);
    end
  end

  for i = 1 : k
    for j = 1 : n
      c_center(i,j) = c_center(i,j) / c_size(i);
    end
  end

  for i = 1 : k
    wss(i) = 0.0;
  end

  for i = 1 : m
    ci = c(i);
    for j = 1 : n
      wss(ci) = wss(ci) + ( a(i,j) - c_center(ci,j) )^2;
    end
  end

  critvl_new = 0.0;
  for i = 1 : k
    critvl_new = critvl_new + wss(i);
  end

  inc = critvl_new - critvl;
%
%  Move object I1 from class C2 to class C1.
%  Move object I2 from class C1 to class C2.
%
  c(i1) = c1;
  c(i2) = c2;

  return
end


text:
% IndexToAssignment Convert index to variable assignment.
%
%   A = IndexToAssignment(I, D) converts an index, I, into the .val vector
%   into an assignment over variables with cardinality D. If I is a vector,
%   then the function produces a matrix of assignments, one assignment
%   per row.
%
%   See also AssignmentToIndex.m and SampleFactors.m
%
% Copyright (C) Daphne Koller, Stanford University, 2012

function A = IndexToAssignment(I, D)

D = D(:)'; % ensure that D is a row vector
A = mod(floor(repmat(I(:) - 1, 1, length(D)) ./ repmat(cumprod([1, D(1:end - 1)]), length(I), 1)), ...
        repmat(D, length(I), 1)) + 1;

end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function [varargout] = funname(varargin)

% PLGNDR associated Legendre function
%
% y = plgndr(n,k,x) computes the values of the associated Legendre functions
% of degree N and order K
%
% implemented as MEX file

% the original implementation was based on "Numerical Recipes in C", version 2.0
% but has been replaced with an equvalent function from GNU Scientific Library

% Copyright (C) 2002, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: plgndr.m 2885 2011-02-16 09:41:58Z roboos $

% compile the missing mex file on the fly
% remember the original working directory
pwdir = pwd;

% determine the name and full path of this function
funname = mfilename('fullpath');
mexsrc  = [funname '.c'];
[mexdir, mexname] = fileparts(funname);

try
  % try to compile the mex file on the fly
  warning('trying to compile MEX file from %s', mexsrc);
  cd(mexdir);
  mex(mexsrc);
  cd(pwdir);
  success = true;

catch
  % compilation failed
  disp(lasterr);
  error('could not locate MEX file for %s', mexname);
  cd(pwdir);
  success = false;
end

if success
  % execute the mex file that was juist created
  funname   = mfilename;
  funhandle = str2func(funname);
  [varargout{1:nargout}] = funhandle(varargin{:});
end

author: PatternRecognition

repo: OpenBMI


text:
%BIPOLAR returns an M-by-3 matrix containing a blue-red colormap, in
%	in which red stands for positive, blue stands for negative,
%	and white stands for 0.
%
%  Usage: cmap = bipolar(M, lo, hi, contrast);  or  cmap = bipolar;
%
%  cmap:  output M-by-3 matrix for BIPOLAR colormap.
%  M:	  number of shades in the colormap. By default, it is the
%	  same length as the current colormap.
%  lo:	  the lowest value to represent.
%  hi:	  the highest value to represent.
%
%  Inspired from the LORETA PASCAL program:
%	http://www.unizh.ch/keyinst/NewLORETA
%
%  jimmy@rotman-baycrest.on.ca
%
%----------------------------------------------------------------
function cmap = bipolar(M, lo, hi, contrast)

   if ~exist('contrast','var')
      contrast = 128;
   end

   if ~exist('lo','var')
      lo = -1;
   end

   if ~exist('hi','var')
      hi = 1;
   end

   if ~exist('M','var')
      cmap = colormap;
      M = size(cmap,1);
   end

   steepness = 10 ^ (1 - (contrast-1)/127);
   pos_infs = 1e-99;
   neg_infs = -1e-99;

   doubleredc = [];
   doublebluec = [];

   if lo >= 0		% all positive

      if lo == 0
         lo = pos_infs;
      end

      for i=linspace(hi/M, hi, M)
         t = exp(log(i/hi)*steepness);
         doubleredc = [doubleredc; [(1-t)+t,(1-t)+0,(1-t)+0]];
      end

      cmap = doubleredc;

   elseif hi <= 0	% all negative

      if hi == 0
         hi = neg_infs;
      end

      for i=linspace(abs(lo)/M, abs(lo), M)
         t = exp(log(i/abs(lo))*steepness);
         doublebluec = [doublebluec; [(1-t)+0,(1-t)+0,(1-t)+t]];
      end

      cmap = flipud(doublebluec);

   else

      if hi > abs(lo)
         maxc = hi;
      else
         maxc = abs(lo);
      end

      for i=linspace(maxc/M, hi, round(M*hi/(hi-lo)))
         t = exp(log(i/maxc)*steepness);
         doubleredc = [doubleredc; [(1-t)+t,(1-t)+0,(1-t)+0]];
      end

      for i=linspace(maxc/M, abs(lo), round(M*abs(lo)/(hi-lo)))
         t = exp(log(i/maxc)*steepness);
         doublebluec = [doublebluec; [(1-t)+0,(1-t)+0,(1-t)+t]];
      end

      cmap = [flipud(doublebluec); doubleredc];

   end

   return;					% bipolar

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Max manipulability index ALONG A LINE.
% Use stomp like to optimize along a surface/line
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function experiment1B
close all;
global robot
global parameters
global hfigures

%STOMP PARAMETERS
%conversion from cost to Prob factor
parameters.lambda = .4;
parameters.lambda_obstacles = .2;
%height of the obstacle
parameters.yo = 2.5;
%cost function starts at this distance
%must be below 0.3 for the 4 DOF robot
parameters.epsilon = 0.2;
%multiply noise by this facto
%parameters.noise_k = 5;
%parameters.noise_sigma_null_space = 0.01;
parameters.alpha=0.02;
parameters.time_step=0.01;

%number of waypoints
parameters.N = 12;
%number of particles
parameters.K = 5;
parameters.n_repeat = 30;
parameters.experiment_name = 'experiment1B_K5_N30.mat';

parameters.animate = 0;

parameters.obstacles = [];

%LINE 1
x1 = -1.5;
y1 = .5; %m
x2 = 0;
y2 = 2; %m
phi = 3*pi/4;
p0 = [x1 y1 0]';
pf = [x2 y2 0]';
T0 = build_T_4dof(p0, phi);
parameters.obstacles{1}.line = [p0 pf];
parameters.obstacles{1}.T0 = T0;

%repeat the experiment E times
random_manips=[];
Gout = [];
for i=1:parameters.n_repeat
    close all
    [pk, final_manip] = stomp_null_space(robot);
    Gout{i}=pk;
    random_manips = [random_manips; final_manip];
    save(parameters.experiment_name)
end

function T = build_T_4dof(p, phi)
T = [cos(phi) -sin(phi) 0 p(1);
     sin(phi) cos(phi) 0 p(2);
     0            0     1  p(3);
     0             0    0   1];

 function T = build_T_sawyer(p, phi)
T = [1  -sin(phi) 0 p(1);
     0  cos(phi) 0 p(2);
     0            0     1  p(3);
     0             0    0   1];

repo: ARTE


text:
function value = angle_rad_nd ( dim_num, v1, v2 )

%*****************************************************************************80
%
%% ANGLE_RAD_ND returns the angle in radians between two rays in ND.
%
%  Discussion:
%
%    This routine always computes the SMALLER of the two angles between
%    two rays.  Thus, if the rays make an (exterior) angle of 1.5 PI,
%    then the (interior) angle of 0.5 PI is reported.
%
%    X dot Y = Norm(X) * Norm(Y) * Cos( Angle(X,Y) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, real V1(DIM_NUM), V2(DIM_NUM), the two rays.
%
%    Output, real ANGLE_RAD_ND, the angle between the rays,
%    in radians.  This value will always be between 0 and PI.
%
  dot = v1(1:dim_num) * v2(1:dim_num)';

  v1norm = sqrt ( sum ( v1(1:dim_num).^2 ) );

  if ( v1norm == 0.0 )
    value = 0.0;
    return
  end

  v2norm = sqrt ( sum ( v2(1:dim_num).^2 ) );

  if ( v2norm == 0.0 )
    value = 0.0;
    return
  end

  value = r8_acos ( dot / ( v1norm * v2norm ) );

  return
end

author: johannesgerer