text:
%% check Clebsch Gordan Tensor

%% the reference

% some arbitrary rotation
g = rotation.byEuler(-72*degree,88*degree,134*degree);

% the rotation matrix
R = matrix(g);

% we want to express the product of two of those rotation matrice

RR_ref = R(:) * R(:).'

%% Express R by D

D1 = WignerD(g,'order',1);
D = D1(:) * D1(:).';

EinsteinSum(U,[1 -1],D1,[-1 -2],conj(U),[2 -2])

%% expansion into Wigner functions of lower order

% zero order component
D0 = WignerD(g,'order',0);
CG0 = ClebschGordanTensor(0);
C0 = D0*EinsteinSum(CG0,[1 3],CG0,[2 4])

% first order component
D1 = WignerD(g,'order',1);
CG1 = ClebschGordanTensor(1);
C1 = EinsteinSum(CG1,[1 3 -1],D1,[-1 -2],CG1,[2 4 -2])

% second order component
D2 = WignerD(g,'order',2);
CG2 = ClebschGordanTensor(2);
C2 = EinsteinSum(CG2,[1 3 -1],D2,[-1 -2],CG2,[2 4 -2])

C = EinsteinSum(C0 + C1 + C2,[-1 -2 -3 -4],U,[1 -1],conj(U),[2 -2],U,[3 -3],conj(U),[4 -4]);
real(reshape(matrix(C),[9 9]))

assert(norm(D - reshape(matrix(C0 + C1 + C2),[9,9]))<=1e-10,'Clebsch Gordan check failed')

%%

% zero order component
D0 = WignerD(g,'order',0);
CG0 = EinsteinSum(ClebschGordanTensor(0),[-1 -2],U,[1 -1],U,[2 -2]);
CG0c = EinsteinSum(ClebschGordanTensor(0),[-1 -2],conj(U),[1 -1],conj(U),[2 -2]);
C0 = D0*EinsteinSum(CG0,[1 3],CG0c,[2 4])

% first order component
D1 = WignerD(g,'order',1);
CG1 = EinsteinSum(ClebschGordanTensor(1),[-1 -2 3],U,[1 -1],U,[2 -2]);
CG1c = EinsteinSum(ClebschGordanTensor(1),[-1 -2 3],conj(U),[1 -1],conj(U),[2 -2]);
C1 = EinsteinSum(CG1,[1 3 -1],D1,[-1 -2],CG1c,[2 4 -2])

% second order component
D2 = WignerD(g,'order',2);
CG2 = EinsteinSum(ClebschGordanTensor(2),[-1 -2 3],U,[1 -1],U,[2 -2]);
CG2c = EinsteinSum(ClebschGordanTensor(2),[-1 -2 3],conj(U),[1 -1],conj(U),[2 -2]);
C2 = EinsteinSum(CG2,[1 3 -1],D2,[-1 -2],CG2c,[2 4 -2])

C = C0 + C1 + C2;
TT_ref ./ real(reshape(matrix(C),[9 9]))

assert(norm(TT_ref - reshape(matrix(C),[9,9]))<=1e-10,'Clebsch Gordan check failed')

author: mtex-toolbox

repo: mtex

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


text:
mkdir image_256
p = dir('./image/');
m1 = [];
for i = 3:numel(p)  %remove . & ..
    str = [ './image/',p(i).name,'/*.jpg'];
    pp = dir(str);
    if(numel(pp)==0)
        str = [ './image/',p(i).name,'/*.bmp'];
        pp = dir(str);
    end
    if(numel(pp)==0)
        str = [ './image/',p(i).name,'/*.png'];
        pp = dir(str);
    end
    for j = 1:numel(pp)
        img_str = [ './image/',p(i).name,'/',pp(j).name];
        wimg_str = [ './image_256/',p(i).name,'/',pp(j).name(1:end-3),'jpg'];
        wdir_str = ['./image_256/',p(i).name];
        mkdir(wdir_str);
        im = imread(img_str);
        sz = size(im);
        if(numel(sz)==2)
           im = repmat(im,1,1,3);
        end
        if(sz(1)<sz(2))
           im = imresize(im,[256, round(sz(2)/sz(1)*256)]);
        else
           im = imresize(im,[round(sz(1)/sz(2)*256),256]);
        end
        imwrite(im,wimg_str);
        m = mean(mean(im,1),2);
        m1 = cat(4,m1,m);
    end
end

mm = mean(m1,4);

author: layumi

repo: Image-Text-Embedding


text:
%%***********************************************************
%% lmiexamp2: generate SDP data for the following LMI problem
%%
%%  max  -Tr(P)
%%  s.t. A1*P + P*A1' + B*diag(d)*B' <= 0
%%       A2*P + P*A2' + B*diag(d)*B' <= 0
%%                                -d <= 0
%%                            sum(d)  = 1
%%***********************************************************
%% Here is an example on how to use this function to
%% find an optimal P.
%%
%% A1 = [-1  0  0; 0 -2  0; 1  1 -1];
%% A2 = A1 + 0.1*A1';
%% B  = [1 2; 3 4; 5 6];
%%
%% [blk,Avec,C,b] = lmiexamp2(A1,A2,B);
%% [obj,X,y,Z] = sqlp(blk,Avec,C,b);
%% n = size(A1,2); N = n*(n+1)/2; dlen = size(B,2);
%% P = smat(blk(1,:),y(1:N));
%% d = y(N+[1:dlen]);
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function [blk,Avec,C,b] = lmiexamp2(A1,A2,B);

%%
   n1 = size(A1,2);  n2 = size(A2,2);
   if (n1 ~= n2); error('lmiexamp2: A1, A2 not compatible'); end;
%%
   n = n1;
   I = speye(n);
   blk{1,1}  = 's'; blk{1,2} = n;
   Avec{1,1} = lmifun(A1,I,B);
   C{1,1} = sparse(n,n);
%%
   blk{2,1}  = 's'; blk{2,2} = n;
   Avec{2,1} = lmifun(A2,I,B);
   C{2,1} = sparse(n,n);
%%
%% and constraints:  -d <= 0
%%
   N = n*(n+1)/2; dlen = size(B,2);
   blk{3,1}  = 'l'; blk{3,2} = dlen;
   Avec{3,1} = [sparse(dlen,N)    -speye(dlen,dlen)];
   C{3,1} = zeros(dlen,1);
%%
%% add in the constraint:  sum(d) = 1
%%
   blk{4,1}  = 'u'; blk{4,2} = 1;
   Avec{4,1} = sparse([zeros(1,N) ones(1,dlen)]);
   C{4,1} = 1;
%%
   blktmp{1,1} = 's'; blktmp{1,2} = n;
   b = [-svec(blktmp,I); zeros(dlen,1)];
%%**********************************************************

author: yu-jiang

repo: radpbook

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


text:
classdef NSGAIIconflict < ALGORITHM
% <many> <real/integer/label/binary/permutation>
% NSGA-II with conflict-based partitioning strategy
% NS     ---  2 --- Number of subspaces
% cycles --- 10 --- Number of cycles

%------------------------------- Reference --------------------------------
% A. L. Jaimes, C. A. Coello Coello, H. Aguirre, and K. Tanaka, Objective
% space partitioning using conflict information for solving many-objective
% problems, Information Sciences, 2014, 268: 305-327.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            [NS,cycles] = Algorithm.ParameterSet(2,10);
            Gc          = ceil(Problem.maxFE/Problem.N/cycles);

            %% Generate random population
            Psi        = {1:Problem.M};
            Population = Problem.Initialization();
            [~,FrontNo,CrowdDis] = EnvironmentalSelection(Population,Problem.N,Psi);

            %% Optimization
            phase = true;
            while Algorithm.NotTerminated(Population)
                MatingPool = TournamentSelection(2,Problem.N,FrontNo,-CrowdDis);
                Offspring  = OperatorGA(Problem,Population(MatingPool));
                [Population,FrontNo,CrowdDis] = EnvironmentalSelection([Population,Offspring],Problem.N,Psi);
                if ~phase && mod(ceil(Problem.FE/Problem.N),Gc)/Gc < 0.3
                    % Change to the approximation phase
                    Psi   = {1:Problem.M};
                    phase = true;
                elseif phase && mod(ceil(Problem.FE/Problem.N),Gc)/Gc >= 0.3
                    % Change to the partitioning phase
                    Psi   = ConflictPartition(Population.objs,NS);
                    phase = false;
                end
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function S2G = delete(S2G,points)
% elilinates points from grid
%
% Input
%  S2G    - @S2Grid
%  points - logical | double | @vector3d
%
% Output
%  S2G    - @S2Grid
%
% See also
% S2Grid/copy

if isa(points,'vector3d'), points = find(S2G,points); end

if isnumeric(points),
  inds = false(length(S2G),1);
  inds(points) = true;
  points = inds;
end

S2G.x(points) = [];
S2G.y(points) = [];
S2G.z(points) = [];

cs = [0 cumsum(GridLength(S2G.rhoGrid))];
% update rho indexing
for ith = 1:GridLength(S2G.thetaGrid)
  irh = points(cs(ith)+1:cs(ith+1));
  S2G.rhoGrid(ith) = delete(S2G.rhoGrid(ith),irh);
end

% update theta indexing
S2G.thetaGrid = delete(S2G.thetaGrid,GridLength(S2G.rhoGrid)==0);
S2G.rhoGrid(GridLength(S2G.rhoGrid) == 0) = [];


text:
%% Lists and Indexing
%
%%
%
% <https://de.mathworks.com/help/matlab/learn_matlab/matrices-and-arrays.html
% Matlab documentation>


text:
function [x,esq,j] = v_kmeanlbg(d,k)
%V_KMEANLBG Vector quantisation using the Linde-Buzo-Gray algorithm [X,ESQ,J]=(D,K)
%
%Inputs:
% D contains data vectors (one per row)
% K is number of centres required
%
%Outputs:
% X is output row vectors (K rows)
% ESQ is mean square error
% J indicates which centre each data vector belongs to
%
%  Implements LBG K-means algorithm:
% Linde, Y., A. Buzo, and R. M. Gray,
% "An Algorithm for vector quantiser design,"
% IEEE Trans Communications, vol. 28, pp.84-95, Jan 1980.

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: v_kmeanlbg.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nc=size(d,2);
[x,esq,j]=v_kmeans(d,1);
m=1;
while m<k
   n=min(m,k-m);
   m=m+n;
   e=1e-4*sqrt(esq)*rand(1,nc);
   [x,esq,j]=v_kmeans(d,m,[x(1:n,:)+e(ones(n,1),:); x(1:n,:)-e(ones(n,1),:); x(n+1:m-n,:)]);
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
%% DEMO_febio_0003_beam_bending
% Below is a demonstration for:
%
% * Building geometry for a beam with hexahedral elements
% * Defining the boundary conditions
% * Coding the febio structure
% * Running the model
% * Importing and visualizing the displacement results

%% Keywords
%
% * febio_spec version 4.0
% * febio, FEBio
% * beam pressure loading
% * surface pressure boundary condition
% * hexahedral elements, hex8
% * beam, rectangular
% * static, solid
% * hyperelastic, Ogden
% * displacement logfile
% * force logfile

%%

clear; close all; clc;

%% Plot settings
fontSize=20;
faceAlpha1=0.8;
markerSize=40;
markerSize2=20;
lineWidth=3;

%% Control parameters

% Path names
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');

% Defining file names
febioFebFileNamePart='tempModel';
febioFebFileName=fullfile(savePath,[febioFebFileNamePart,'.feb']); %FEB file name
febioLogFileName=[febioFebFileNamePart,'.txt']; %FEBio log file name
febioLogFileName_disp=[febioFebFileNamePart,'_disp_out.txt']; %Log file name for exporting displacement
febioLogFileName_force=[febioFebFileNamePart,'_force_out.txt']; %Log file name for exporting force

%Specifying dimensions and number of elements
beamWidth=10;
sampleWidth=beamWidth; %Width
sampleThickness=4*beamWidth; %Thickness
sampleHeight=beamWidth; %Height
pointSpacings=2*ones(1,3); %Desired point spacing between nodes
numElementsWidth=round(sampleWidth/pointSpacings(1)); %Number of elemens in dir 1
numElementsThickness=round(sampleThickness/pointSpacings(2)); %Number of elemens in dir 2
numElementsHeight=round(sampleHeight/pointSpacings(3)); %Number of elemens in dir 3

%Define applied force
appliedPressure=1e-5;

%Material parameter set
c1=1e-3; %Shear-modulus-like parameter
m1=8; %Material parameter setting degree of non-linearity
k_factor=1e2; %Bulk modulus factor
k=c1*k_factor; %Bulk modulus

% FEA control settings
numTimeSteps=10; %Number of time steps desired
max_refs=25; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=6; %Optimum number of iterations
max_retries=5; %Maximum number of retires
dtmin=(1/numTimeSteps)/100; %Minimum time step size
dtmax=1/numTimeSteps; %Maximum time step size

runMode='external';% 'internal' or 'external'

%% Creating model geometry and mesh
% A box is created with tri-linear hexahedral (hex8) elements using the
% |hexMeshBox| function. The function offers the boundary faces with
% seperate labels for the top, bottom, left, right, front, and back sides.
% As such these can be used to define boundary conditions on the exterior.

% Create a box with hexahedral elements
beamDimensions=[sampleWidth sampleThickness sampleHeight]; %Dimensions
beamElementNumbers=[numElementsWidth numElementsThickness numElementsHeight]; %Number of elements
outputStructType=2; %A structure compatible with mesh view
[meshStruct]=hexMeshBox(beamDimensions,beamElementNumbers,outputStructType);

%Access elements, nodes, and faces from the structure
E=meshStruct.elements; %The elements
V=meshStruct.nodes; %The nodes (vertices)
Fb=meshStruct.facesBoundary; %The boundary faces
Cb=meshStruct.boundaryMarker; %The "colors" or labels for the boundary faces
elementMaterialIndices=ones(size(E,1),1); %Element material indices

%%
% Plotting model boundary surfaces and a cut view

hFig=cFigure;

subplot(1,2,1); hold on;
title('Model boundary surfaces and labels','FontSize',fontSize);
gpatch(Fb,V,Cb,'k',faceAlpha1);
colormap(gjet(6)); icolorbar;
axisGeom(gca,fontSize);

hs=subplot(1,2,2); hold on;
title('Cut view of solid mesh','FontSize',fontSize);
optionStruct.hFig=[hFig hs];
meshView(meshStruct,optionStruct);
axisGeom(gca,fontSize);

drawnow;

%% Defining the boundary conditions
% The visualization of the model boundary shows colors for each side of the
% cube. These labels can be used to define boundary conditions.

%Define supported node set
bcSupportList=unique(Fb(Cb==4,:)); %Node set part of selected face

%Loaded surface
F_pressure=Fb(Cb==6,:);

%%
% Visualizing boundary conditions. Markers plotted on the semi-transparent
% model denote the nodes in the various boundary condition lists.

hf=cFigure;
title('Boundary conditions','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;

gpatch(Fb,V,'kw','k',0.5);

hl(1)=plotV(V(bcSupportList,:),'k.','MarkerSize',markerSize);
hl(2)=gpatch(F_pressure,V,'r','k',1);

legend(hl,{'BC support','BC prescribe'});

axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Defining the FEBio input structure
% See also |febioStructTemplate| and |febioStruct2xml| and the FEBio user
% manual.

%Get a template with default settings
[febio_spec]=febioStructTemplate;

%febio_spec version
febio_spec.ATTR.version='4.0';

%Module section
febio_spec.Module.ATTR.type='solid';

%Control section
febio_spec.Control.analysis='STATIC';
febio_spec.Control.time_steps=numTimeSteps;
febio_spec.Control.step_size=1/numTimeSteps;
febio_spec.Control.solver.max_refs=max_refs;
febio_spec.Control.solver.qn_method.max_ups=max_ups;
febio_spec.Control.time_stepper.dtmin=dtmin;
febio_spec.Control.time_stepper.dtmax=dtmax;
febio_spec.Control.time_stepper.max_retries=max_retries;
febio_spec.Control.time_stepper.opt_iter=opt_iter;

%Material section
materialName1='Material1';
febio_spec.Material.material{1}.ATTR.name=materialName1;
febio_spec.Material.material{1}.ATTR.type='Ogden';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.c1=c1;
febio_spec.Material.material{1}.m1=m1;
febio_spec.Material.material{1}.c2=c1;
febio_spec.Material.material{1}.m2=-m1;
febio_spec.Material.material{1}.k=k;

%Mesh section
% -> Nodes
febio_spec.Mesh.Nodes{1}.ATTR.name='nodeSet_all'; %The node set name
febio_spec.Mesh.Nodes{1}.node.ATTR.id=(1:size(V,1))'; %The node id's
febio_spec.Mesh.Nodes{1}.node.VAL=V; %The nodel coordinates

% -> Elements
partName1='Part1';
febio_spec.Mesh.Elements{1}.ATTR.name=partName1; %Name of this part
febio_spec.Mesh.Elements{1}.ATTR.type='hex8'; %Element type
febio_spec.Mesh.Elements{1}.elem.ATTR.id=(1:1:size(E,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=E; %The element matrix

% -> Surfaces
surfaceName1='LoadedSurface';
febio_spec.Mesh.Surface{1}.ATTR.name=surfaceName1;
febio_spec.Mesh.Surface{1}.quad4.ATTR.id=(1:1:size(F_pressure,1))';
febio_spec.Mesh.Surface{1}.quad4.VAL=F_pressure;

% -> NodeSets
nodeSetName1='bcSupportList';
febio_spec.Mesh.NodeSet{1}.ATTR.name=nodeSetName1;
febio_spec.Mesh.NodeSet{1}.VAL=bcSupportList(:)';

%MeshDomains section
febio_spec.MeshDomains.SolidDomain.ATTR.name=partName1;
febio_spec.MeshDomains.SolidDomain.ATTR.mat=materialName1;

%Boundary condition section
% -> Fix boundary conditions
febio_spec.Boundary.bc{1}.ATTR.name='zero_displacement_xyz';
febio_spec.Boundary.bc{1}.ATTR.type='zero displacement';
febio_spec.Boundary.bc{1}.ATTR.node_set=nodeSetName1;
febio_spec.Boundary.bc{1}.x_dof=1;
febio_spec.Boundary.bc{1}.y_dof=1;
febio_spec.Boundary.bc{1}.z_dof=1;

%Loads section
% -> Surface load
febio_spec.Loads.surface_load{1}.ATTR.type='pressure';
febio_spec.Loads.surface_load{1}.ATTR.surface=surfaceName1;
febio_spec.Loads.surface_load{1}.pressure.ATTR.lc=1;
febio_spec.Loads.surface_load{1}.pressure.VAL=appliedPressure;
febio_spec.Loads.surface_load{1}.symmetric_stiffness=1;

%LoadData section
% -> load_controller
febio_spec.LoadData.load_controller{1}.ATTR.id=1;
febio_spec.LoadData.load_controller{1}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{1}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{1}.points.point.VAL=[0 0; 1 1];

%Output section
% -> log file
febio_spec.Output.logfile.ATTR.file=febioLogFileName;
febio_spec.Output.logfile.node_data{1}.ATTR.file=febioLogFileName_disp;
febio_spec.Output.logfile.node_data{1}.ATTR.data='ux;uy;uz';
febio_spec.Output.logfile.node_data{1}.ATTR.delim=',';

febio_spec.Output.logfile.node_data{2}.ATTR.file=febioLogFileName_force;
febio_spec.Output.logfile.node_data{2}.ATTR.data='Rx;Ry;Rz';
febio_spec.Output.logfile.node_data{2}.ATTR.delim=',';

%% Quick viewing of the FEBio input file structure
% The |febView| function can be used to view the xml structure in a MATLAB
% figure window.

%%
% |febView(febio_spec); %Viewing the febio file|

%% Exporting the FEBio input file
% Exporting the febio_spec structure to an FEBio input file is done using
% the |febioStruct2xml| function.

febioStruct2xml(febio_spec,febioFebFileName); %Exporting to file and domNode

%% Running the FEBio analysis
% To run the analysis defined by the created FEBio input file the
% |runMonitorFEBio| function is used. The input for this function is a
% structure defining job settings e.g. the FEBio input file name. The
% optional output runFlag informs the user if the analysis was run
% succesfully.

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.disp_log_on=1; %Display convergence information in the command window
febioAnalysis.runMode=runMode;

[runFlag]=runMonitorFEBio(febioAnalysis);%START FEBio NOW!!!!!!!!

%% Import FEBio results

if runFlag==1 %i.e. a succesful run

    %%
    % Importing nodal displacements from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_disp),1,1);

    %Access data
    N_disp_mat=dataStruct.data; %Displacement
    timeVec=dataStruct.time; %Time

    %Create deformed coordinate set
    V_DEF=N_disp_mat+repmat(V,[1 1 size(N_disp_mat,3)]);

    %%
    % Plotting the simulated results using |anim8| to visualize and animate
    % deformations

    DN_magnitude=sqrt(sum(N_disp_mat(:,:,end).^2,2)); %Current displacement magnitude

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    title('Displacement magnitude [mm]','Interpreter','Latex')
    hp=gpatch(Fb,V_DEF(:,:,end),DN_magnitude,'k',1); %Add graphics object to animate
    hp.Marker='.';
    hp.MarkerSize=markerSize2;
    hp.FaceColor='interp';
    gpatch(Fb,V,0.5*ones(1,3),'k',0.25); %A static graphics object

    axisGeom(gca,fontSize);
    colormap(gjet(250)); colorbar;
    caxis([0 max(DN_magnitude)]);
    axis(axisLim(V_DEF)); %Set axis limits statically
        camlight headlight;

    % Set up animation features
    animStruct.Time=timeVec; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments
        DN_magnitude=sqrt(sum(N_disp_mat(:,:,qt).^2,2)); %Current displacement magnitude

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp hp]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData'}; %Properties of objects to animate
        animStruct.Set{qt}={V_DEF(:,:,qt),DN_magnitude}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature
    drawnow;

end

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
classdef MimEditBlue < MimGuiPlugin
    % MimEditBlue. Gui Plugin for setting paint colour
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimEditBlue is a Gui Plugin for the MIM Toolkit.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Blue'
        SelectedText = 'Blue'
        ToolTip = 'Edit with blue label'
        Category = 'Segmentation label'
        Visibility = 'Dataset'
        Mode = 'Edit'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 5
        ButtonHeight = 1
        Icon = 'paint.png'
        IconColour = GemMarkerPoint.DefaultColours{1}
        Location = 31
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.PaintBrushColour = 1;
        end

        function enabled = IsEnabled(gui_app)
             enabled = gui_app.IsDatasetLoaded && gui_app.ImagePanel.OverlayImage.ImageExists && ...
                isequal(gui_app.ImagePanel.SelectedControl, 'Paint');
        end

        function is_selected = IsSelected(gui_app)
            is_selected = gui_app.ImagePanel.PaintBrushColour == 1;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function c = appendevent(a, b)

% APPENDEVENT

% Copyright (C) 2008, 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: appendevent.m 2885 2011-02-16 09:41:58Z roboos $

if isempty(a)
  c = b(:);
elseif isempty(b)
  c = a(:);
else
  c = a(:);
  for i=1:numel(b)
    c(end+1).type     = b(i).type;
    c(end  ).value    = b(i).value;
    c(end  ).sample   = b(i).sample;
    if isfield(b, 'timestamp')
      c(end  ).timestamp = b(i).timestamp; % optional
    end
    if isfield(b, 'offset')
      c(end  ).offset    = b(i).offset;    % optional
    end
    if isfield(b, 'duration')
      c(end  ).duration  = b(i).duration;  % optional
    end
  end
end

author: PatternRecognition

repo: OpenBMI


text:
classdef HButtonBox < uix.HButtonBox
    %uiextras.HButtonBox  Arrange buttons horizontally in a single row
    %
    %   obj = uiextras.HButtonBox() is a type of HBox specialised for
    %   arranging a row of buttons, check-boxes or similar graphical
    %   elements. All buttons are given equal size and by default are
    %   centered in the drawing area. The justification can be changed as
    %   required.
    %
    %   obj = uiextras.HButtonBox(param,value,...) also sets one or more
    %   parameter values.
    %
    %   See the documentation for more detail and the list of properties.
    %
    %   Examples:
    %   >> f = figure();
    %   >> b = uiextras.HButtonBox( 'Parent', f );
    %   >> uicontrol( 'Parent', b, 'String', 'One' );
    %   >> uicontrol( 'Parent', b, 'String', 'Two' );
    %   >> uicontrol( 'Parent', b, 'String', 'Three' );
    %   >> set( b, 'ButtonSize', [130 35], 'Spacing', 5 );
    %
    %   See also: uiextras.VButtonBox
    %             uiextras.HBox

    %  Copyright 2009-2014 The MathWorks, Inc.
    %  $Revision: 1077 $ $Date: 2015-03-19 16:44:14 +0000 (Thu, 19 Mar 2015) $

    methods

        function obj = HButtonBox( varargin )
            %uiextras.HButtonBox  Create a new horizontal button box

            % Call uix constructor
            obj@uix.HButtonBox( varargin{:} )

            % Auto-parent
            if ~ismember( 'Parent', varargin(1:2:end) )
                obj.Parent = gcf();
            end

        end % constructor

    end % structors

end % classdef

author: goGPS-Project


text:
function [W,dist_H_V] = shepard(V,C,P,E,CE,p)
  % SHEPARD Compute shepard weights for a list of vertices, given a list of
  % samples and optionally a denominator power value.
  %
  % W = shepard(V,C)
  %
  % Inputs:
  %  V  list of vertex positions
  %  C  list of control vertices
  %  P  list of indices into C for point controls, { 1:size(C,1) }
  %  E  list of bones, pairs of indices into C, connecting control vertices,
  %    { [] }
  %  CE  list of "cage edges", pairs of indices into ***P***, connecting
  %    control ***points***. A "cage edge" just tells point boundary conditions
  %    to vary linearly along straight lines between the end points, and to be
  %    zero for all other handles. { [] }
  %  p  (optional) power for denominator, scalar or list same size as C {2}
  %
  % Outputs:
  %  W  weights, # vertices by # handles matrix of weights
  %
  % Example:
  %  % "Shape-aware" shepard weights
  %  D = permute(sum((repmat(V,[1,1,size(C,1)]) - ...
  %    permute(repmat(C,[1,1,n]),[3,2,1])).^2,2),[1,3,2]);
  %  [minD,b] = min(D);
  %  snap_C = V(b,:);
  %  [BV,ev,ed] = biharmonic_embedding(V,F,6);
  %  [W] = shepard(BV,BV(b,:));
  %
  %

  % check if either are 3D but really all z's are 0
  V_flat = size(V,2) == 3 && (sqrt(sum(V(:,3).^2)) < 1e-10);
  C_flat = size(C,2) == 3 && (sqrt(sum(C(:,3).^2)) < 1e-10);
  % is both are essentially 2D then ignore z-coords
  if((size(C,2) == 2 || C_flat) && (size(V,2) == 2 || V_flat))
    % ignore z coordinate
    V = V(:,1:2);
    C = C(:,1:2);
  end

  assert(size(C,2) == size(V,2));
  dim = size(C,2);

  % default p value
  if(~exist('p','var'))
    p = 2;
  end

  if(~exist('P','var'))
    P = 1:size(C,1);
  end

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

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

  assert(prod(size(P)) == max(size(P)));
  assert(isempty(E) || size(E,2) == 2);
  assert(isempty(CE) || size(CE,2) == 2);
  assert(isempty(E) || max(E(:))<=size(C,1));
  assert(isempty(CE) || max(CE(:))<=numel(P));
  assert(isempty(P) || max(P(:))<=size(C,1));
  assert(isempty(E) || min(E(:)) >= 1);
  assert(isempty(CE) || min(CE(:)) >= 1);
  assert(isempty(P) || min(P(:)) >= 1);

  % number of point controls
  np = numel(P);
  nb = size(E,1);
  % number of domain vertices
  n = size(V,1);

  % if p is a scalar convert it now to a list the same size as C
  if(prod(size(p)) == 1)
    p = repmat(p,np+nb,1);
  elseif ~isempty(CE)
    if 0 ~= std(p(CE(:)))
      error('Variable powers not supported with cage edges');
    end
  end

  dist_P_V = zeros(np,n);
  if(np > 0)
    % vectors from V to every P, where PmV(i,j,:) is the vector from domain
    % vertex j to handle i
    PmV = ...
      permute( ...
        permute(repmat(C(P,:),[1,1,n]),[3,2,1]) - ...
        repmat(V,[1,1,np]),[3,1,2]);
    % distance from V to every P, where dist_P_V(i,j) is the distance from domain
    % vertex j to point handle i
    dist_P_V = sqrt(sum(PmV.^2,3));
    % distance from each corner in P to the next corner so that edge_length(i)
  end

  dist_B_V = zeros(nb,n);
  if(nb > 0)
    % loop over bones
    for( ii = 1:nb )
      % project to line
      [t,sqr_d] = project_to_lines(V,C(E(ii,1),:),C(E(ii,2),:));
      d = sqrt(sqr_d);
      % if projected point falls outside of line segment take distance to
      % closest end point
      d(t<0) = sqrt(sum((V(t<0,:)-repmat(C(E(ii,1),:),size(V(t<0,:),1),1)).^2,2));
      d(t>1) = sqrt(sum((V(t>1,:)-repmat(C(E(ii,2),:),size(V(t>1,:),1),1)).^2,2));
      dist_B_V(ii,:) = d';
    end
  end

  dist_H_V = [dist_P_V ; dist_B_V];
  % power of each control point seen by each vertex in domain
  pp = repmat(p,1,n);
  W = 1.0./((dist_H_V).^pp);

  if(size(CE,1) > 0)
    % zero out previous entries in points with incident cage edges
    W(CE(:),:) = 0;
    % compute projection of each point to each line segment
    [T,sqrD] = project_to_lines(V,C(P(CE(:,1)),:),C(P(CE(:,2)),:));
    % compute weights for each cage edge as if it were a bone
    pce = (p(CE(:,1))+p(CE(:,2)))/2;
    [~,dist_CE_V] = shepard(V,C,[],P(CE),[],pce);
    % compute "bone weights" for each cage edge as they were before
    % normalization
    ppce = repmat(pce,1,n);
    WCE = 1.0./((dist_CE_V).^ppce)';
    % clamp distances
    T(T<0) = 0;
    T(T>1) = 1;
    % multiply "bone weight" by projection parameter and sum up
    II = repmat(1:n,1,2*size(CE,1))';
    JJ = reshape(repmat(CE(:)',n,1),[],1);
    VV = [WCE(:).*(1-T(:)); WCE(:).*T(:)];
    W = W+sparse(II,JJ,VV,n,np+nb)';
    % cage edges on-samples
    on_sample = dist_CE_V < eps;
    [I,J] = find(on_sample);
    W(:,any(on_sample)) = 0;
    % grab T transpose
    TT = T';
    W(sub2ind(size(W),CE(I,1),J)) = 1-TT(on_sample);
    W(sub2ind(size(W),CE(I,2),J)) = TT(on_sample);
  end

  % Handle degenrate case that a control point is on a mesh vertex
  % snap vertices close to corners
  on_sample = dist_H_V < eps;
  W(:,any(on_sample,1)) = 0;
  W(on_sample) = 1;

  % normalize W
  W = W./repmat(sum(W,1),np+nb,1);

  % we've made W transpose
  W = W';

end

author: alecjacobson

repo: gptoolbox


text:
function c8_uniform_01_test ( )

%*****************************************************************************80
%
%% C8_UNIFORM_01_TEST tests C8_UNIFORM_01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 June 2006
%
%  Author:
%
%    John Burkardt
%
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8_UNIFORM_01_TEST\n' );
  fprintf ( 1, '  C8_UNIFORM_01 computes pseudorandom complex values\n' );
  fprintf ( 1, '  in the unit circle.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The initial seed is %d\n', seed );
  fprintf ( 1, '\n' );

  for i = 1 : 10
    [ x, seed ] = c8_uniform_01 ( seed );
    fprintf ( 1, '  %6d  ( %f, %f )\n', i, real ( x ), imag ( x ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Example A.1: Linear Advection in a Periodic Domain
%
% The linear advection problem with periodic boundary conditions
%
% $$u_t + u_x = 0, \qquad u(x,0)=u_0(x), \qquad u(0,t)=u(1,t)$$
%
% is well suited for studying error mechanisms in numerical schemes for
% hyperbolic conservation laws. In particular, to study dissipative and
% oscillatory errors, we will use initial data consisting of a combination
% of a smooth, squared cosine wave and double step function.

%%
% We consider two first-order methods (Lax-Friedrichs and the upwind
% method) and two second-order methods (Lax-Wendroff and MacCormack)
T  = 10;
N  = 100; h=1/N; x=h*(1:N);
u0 = (abs(x-.25)<=0.15).*(cos(pi*(10/3)*(x-0.25))).^2 + ...
   1.0*(abs(x-0.75)<0.15);
xx  = linspace(0,1,1001);
uf  = (abs(xx-.25)<=0.15).*(cos(pi*(10/3)*(xx-0.25))).^2 + ...
   1.0*(abs(xx-0.75)<0.15);

method = {'LxF', 'upwind', 'LxW', 'McC'};
name   = {'Lax-Friedrichs', 'Upwind', 'Lax-Wendroff', 'MacCormack'};
for i=1:4
   u=finvol('linearfunc',u0,T,h,1,'periodic',method{i},.9);
   subplot(2,2,i);
   plot(xx,uf,'-',x,u,'o','MarkerSize',2); axis([0 1 -0.2 1.3]);
   title(name{i});
end
%%
% We see that the two first-order schemes smear both the smooth part and
% the discontinuous path of the advected profile and that the
% Lax-Friedrichs method is more diffusive than the upwind method. The
% second-order schemes, on the other hand, preserve the smooth profile
% quite accurately, but introduce spurious oscillations around the two
% discontinuities.

repo: Aero-matlab


text:
function setLineProperty(h)

set(h(2),'LineWidth',2)
set(h(3),'LineStyle',':')

author: s-kajita

repo: IntroductionToHumanoidRobotics


text:
function SO3F = log(SO3F)
% overloads |log(SO3F)|
%
% Syntax
%   SO3F = log(SO3F)
%
% Input
%  SO3F - @SO3Fun
%
% Output
%  SO3F - @SO3Fun
%

SO3F = SO3FunHandle(@(rot) log(SO3F.eval(rot)),SO3F.SRight,SO3F.SLeft);

end


text:
% The function gamer creates the objective function and the constraints for
% the optimization problem to supply it to fmincon.

function [x,fval,exitflag,output] = gamer(n,Us,p,I,s,ub,lb,x0,Aeq,beq,pay,U)

    function F = myfun(x)
        Funct = 0;
        prod = 1;
        for i = 1 : n
            Funct = Funct + x(s+i);
        end
        for i = 1 : p
            for j = 1 : n
                prod = prod * x(I(i,j));
            end
            Funct = Funct - Us(i) * prod;
            prod = 1;
        end
        F = Funct;
    end

    function [c ceq] = confun(x)
        C = zeros(s,1);
        for i = 1 : s
            C(i) = -x(pay(i));
            for t = 1 : n
                add = 0;
                for j = 1 : p
                    prd = 1;
                    for k = 1 : n
                        if i == I(j,k)
                            prd = prd * U(j,k);
                        else
                            prd = prd * x(I(j,k));
                        end
                    end
                    if I(j,t) ~= i
                        prd = 0;
                    end
                    add = add + prd;
                end
                C(i) = add + C(i);
            end
        end
        c = C;
        ceq = [];
    end

options = optimset('Display','off');
warning 'off' 'all';
[x,fval,exitflag,output] = fmincon(@myfun,x0,[],[],Aeq,beq,lb,ub,@confun,options);

end

author: Sable

repo: mcbench-benchmarks


text:
function status = svm_learn(options, examples, model)
% SVM_LEARN - Interface to SVM light, learning module
%
%   STATUS = SVM_LEARN(OPTIONS, EXAMPLES, MODEL)
%   Call the training program 'svm_learn' of the SVM light
%   package.
%   OPTIONS must be a structure generated by SVMLOPT.
%   EXAMPLES is the name of the file containing the training examples
%   (use SVMLWRITE to convert a Matlab matrix to the appropriate format).
%   MODEL is the name of the file holding the trained Support Vector
%   Machine.
%   If 'svm_learn' is not on the path, OPTIONS must contain a field
%   'ExecPath' with the path of the executable.
%   STATUS is the error code returned by SVM light (0 if everything went
%   fine)
%
%   See also SVMLOPT, SVMLWRITE, SVM_CLASSIFY, SVMLREAD
%

%
% Copyright (c) by Anton Schwaighofer (2001)
% $Revision: 1.6 $ $Date: 2002/08/09 20:24:03 $
% mailto:anton.schwaighofer@gmx.net
%
% This program is released unter the GNU General Public License.
%

error(nargchk(3, 3, nargin));

% check parameter consistency for kernels
if ~isempty(options.Kernel),
  if (options.Kernel~=0) & isempty(options.KernelParam),
    error(sprintf('The chosen Kernel = %i requires parameters, but none are given in KernelParam', ...
                  options.Kernel));
  end
  parlen = length(options.KernelParam);
  isString = isa(options.KernelParam, 'char');
  switch options.Kernel
    case {1, 2}
      if isString | (parlen~=1),
        error(sprintf('The chosen Kernel = %i requires a scalar parameter', ...
                      options.Kernel));
      end
    case 3
       if isString | (parlen~=2),
        error(sprintf('The chosen Kernel = %i requires 2 scalar parameters', ...
                      options.Kernel));
      end
    case 4,
      if ~isString,
        error(sprintf('The chosen Kernel = %i requires a string parameter', ...
                      options.Kernel));
      end
  end
end

if ~isempty(options.NewVariables),
  if isempty(options.MaximumQP),
    maxval = 10;
  else
    maxval = options.MaximumQP;
  end
  if options.NewVariables>maxval,
    error('Option ''NewVariables'' must be smaller than 10 resp. value of MaximumQP');
  end
end

Names = fieldnames(options);
[m,n] = size(Names);

s = '';
for i = 1:m,
  field = Names{i,:};
  value = getfield(options, field);
  switch field,
    case 'Verbosity'
      s = stroption(s, '-v %i', value);
    case 'Regression'
      if ~isempty(value),
        if value==0,
          s = [s ' -z c'];
        else
          s = [s ' -z r'];
        end
      end
    case 'C'
      s = stroption(s, '-c %.10g', value);
    case 'TubeWidth'
      s = stroption(s, '-w %.10g', value);
    case 'CostFactor'
      s = stroption(s, '-j %.10g', value);
    case 'Biased'
      s = stroption(s, '-b %i', value);
    case 'RemoveIncons'
      s = stroption(s, '-i %i', value);
    case 'ComputeLOO'
      s = stroption(s, '-x %i', value);
    case 'XialphaRho'
      s = stroption(s, '-o %.10g', value);
    case 'XialphaDepth'
      s = stroption(s, '-k %.10g', value);
    case 'TransPosFrac'
      s = stroption(s, '-p %.10g', value);
    case 'Kernel'
      s = stroption(s, '-t %i', value);
    case 'KernelParam'
      if ~isempty(value),
        switch options.Kernel
          case 0
          case 1
            s = stroption(s, '-d %.10g', value(1));
          case 2
            s = stroption(s, '-g %.10g', value(1));
          case 3
            s = stroption(s, '-s %.10g -r %.10g', value(1), value(2));
          case 4
            s = stroption(s, '-u "%s"', value);
        end
      end
    case 'MaximumQP'
      s = stroption(s, '-q %i', value);
    case 'NewVariables'
      s = stroption(s, '-n %i', value);
    case 'CacheSize'
      s = stroption(s, '-m %i', value);
    case 'EpsTermin'
      s = stroption(s, '-e %.10g', value);
    case 'ShrinkIter'
      s = stroption(s, '-h %i', value);
    case 'ShrinkCheck'
      s = stroption(s, '-f %i', value);
    case 'TransLabelFile'
      s = stroption(s, '-l %s', value);
    case 'AlphaFile'
      s = stroption(s, '-a %s', value);
  end
end

evalstr = [fullfile(options.ExecPath, 'svm_learn') s ' ' ...
           examples ' ' model];
% fprintf('\nCalling SVMlight:\n%s\n\n', evalstr);
if isunix,
  status = unix(evalstr);
else
  status = dos(evalstr);
end

function s = stroption(s, formatstr, value, varargin)
% STROPTION - Add a new option to string
%

if ~isempty(value),
  s = [s ' ' sprintf(formatstr, value, varargin{:})];
end

author: viggin

repo: domain-adaptation-toolbox


text:
function [U,T,mu] = pcasecon(X,k)
% Input:
% X : m x n matrix
% Each column of X is a feature vector
%
% Output:
% X = U*T approximately (up to k)
%
% Description:
% Principal Component Analysis (PCA) while trying to conserve memory and be faster
% Requires that k <= min(m,n) where [m,n] = size(X)
% This function is useful if k is much smaller than m and n
% or if X is sparse (see doc eigs)
%
% Vipin Vijayan (2014)

mu = mean(X,2);
X = bsxfun(@minus,X,mu);

[m,n] = size(X);
assert(k <= m && k <= n, 'k needs to be smaller than size(X,1) and size(X,2)');

if  m <= n
    C = X*X';
    [U,~] = eigs(C,k);
    clear C;

    T = U'*X;
else
    C = X'*X;
    [V,D] = eigs(C,k);
    clear C;

    U = X*V; % convert evecs from X'*X to X*X'. the evals are the same.
    %s = sqrt(sum(U.^2,1))';
    s = sqrt(abs(diag(D)));
    U = bsxfun(@(x,c)x./c, U, s');
    %S = diag(s);
    %T = S*V';
    T = bsxfun(@(c,vt)c.*vt,s,V');
end

author: andrewssobral

repo: lrslibrary


text:
function MRSIStruct = setDimsAllZero(MRSIStruct)
    dimNames = fieldnames(MRSIStruct.dims);
    for iName = 1:length(dimNames)
        MRSIStruct = setDimension(MRSIStruct, dimNames{iName}, 0);
    end
end

author: CIC-methods

repo: FID-A

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


text:
function results = searchModel(model,searchTerm,varargin)
% Search for the specified term in the given model.
% The performed search is fuzzy if similarity is lower than 1.
%
% USAGE:
%    results = searchModel(model,searchTerm,...)
%
% INPUTS:
%    model:         The model to search in
%    searchTerm:    The term to search for
%    varargin:      Additional parameters as parameter/value pairs or value struct.
%                   Available parameters are:
%                    * printLevel - Printout switch, 0 silent, 1 print results. (Default = 1) .
%                    * similarity - Minimum similarity (as provided in calcSim) for matches.
%
% OUTPUTS:
%    results:       The results struct build as follows:
%
%                     * .field - a struct array containing of the basic fields (e.g. rxns, mets etc) of the model in which matches were found with each field indicating the matches as detailed below.
%                     * .field.id - The id of the matching element
%                     * .field.matches - information about the matches of the respective ID.
%                     * .field.matches.value - the value that matched.
%                     * .field.matches.source - the field that contained the matching value.
%
% .. Author: - Thomas Pfau, June 2018

parser = inputParser();
parser.addParameter('printLevel',1,@isnumeric);
parser.addParameter('similarity',0.8,@(x) isnumeric(x) && x <= 1 && x >= 0);
parser.parse(varargin{:});
similarity = parser.Results.similarity;
printLevel = parser.Results.printLevel;

% first collect some potential fields, which are known.
baseFields = {'rxns','mets','vars','ctrs','genes','comps','prots'};
baseFieldNames = {'reactions', 'metabolites', 'variables', 'constraints', 'genes', 'compartments', 'proteins'};

%filter the base fields:
fieldsPresentInModel = ismember(baseFields,fieldnames(model));
baseFields = baseFields(fieldsPresentInModel);
baseFieldNames = baseFieldNames(fieldsPresentInModel);

%All Base fields do have an associated "Names" field.
nameFields = regexprep(baseFields,'s$','Names');
knownFields = union(baseFields,nameFields);
dbFields = getDefinedFieldProperties('DataBaseFields',true);
knownFields = union(knownFields,dbFields(:,3));
% and get the annotations which can also be looked up.
annotationQualifiers = getBioQualifiers();
results = struct();
% now, loop over all basic fields
displayedFields = 0;
for field = 1:numel(baseFields)
    cField = baseFields{field};
    % get the model fields associated with this type.
    modelFields = getModelFieldsForType(model,cField);
    resultList = cell(numel(model.(cField)),numel(modelFields));
    similarities = zeros(numel(model.(cField)),1);
    annotationsFields = {};
    for aqual = 1:numel(annotationQualifiers)
       cAnnotType = regexprep(cField,'s$',annotationQualifiers{aqual});
       annotationsFields = union(annotationsFields,modelFields(cellfun(@(x) strncmp(x,cAnnotType,length(cAnnotType)),modelFields)));
    end
    % except for xyzNames there are few other fields which contain sensibly
    % searchable information (subSystems is one example).
    for modelField = 1:numel(modelFields)
        cModelField = modelFields{modelField};
        if strcmp(cModelField,'subSystems')

            if iscell(model.subSystems{1})
                % subSystems is special, as it contains cell arrays.
                fieldToUse = cellfun(@(x) strjoin(x,';'),model.subSystems,'Uniform',0);
                isAnnotation = true;
            else
                fieldToUse = model.subSystems;
                isAnnotation = false;
            end

        elseif any(strcmp(cModelField,knownFields))
            % iIf its a known field, than it is a cell array of strings and
            % we will use it accordingly
            isAnnotation = false;
            if strcmp(cModelField,'mets')
                if isempty(regexp(searchTerm,'\[[^\[]\]$'))
                    fieldToUse = regexprep(model.mets,'\[[^\[]\]$','');
                else
                    fieldToUse = model.mets;
                end
            else
                fieldToUse = model.(cModelField);
            end
        else
            if any(strcmp(cModelField,annotationsFields))
                % if its an annotation field, we look into it.
                fieldToUse = model.(cModelField);
                isAnnotation = true;
            else
                % field does nto match anything. don't search in it.
                continue
            end
        end
        % find matches
        [matchingIDs,positions,csims] = findMatchingFieldEntries(fieldToUse,searchTerm,isAnnotation,similarity);
        if ~isempty(positions)
            resultList(positions,modelField) = matchingIDs;
            similarities(positions) = max(similarities(positions),csims);
        end
    end
    if any(any(~cellfun(@isempty, resultList)))
        if printLevel > 0 && displayedFields > 0
            fprintf('\n-----------------------------------------------------\n\n');
        end
        displayedFields = displayedFields + 1;
        if printLevel > 0
            fprintf('The following %s have matching properties:\n\n',baseFieldNames{field});
        end
        % get the base field
        results.(cField) = struct();
        relRows = ~all(cellfun(@isempty,resultList),2);
        results.(cField).id = 'start';
        results.(cField).matches = struct();
        results.(cField)(sum(relRows)).id = 'end';
        % and the relevant results for that base field
        relResults = resultList(relRows,:);
        relResultPos = find(relRows);
        % order according to highest similarity
        [~,simorder] = sort(similarities(relRows),'descend');
        for cResults = 1:size(relResults,1)
            results.(cField)(cResults).id = model.(cField){relResultPos(simorder(cResults))};
            results.(cField)(cResults).matches = struct();
            % init struct with size
            results.(cField)(cResults).matches.source = '';
            results.(cField)(cResults).matches.value = '';
            % get the field which have matching entries.
            resultEntries = find(~cellfun(@isempty, relResults(simorder(cResults),:)));
            results.(cField)(cResults).matches(numel(resultEntries)).source = '';
            for cResult = 1:numel(resultEntries)
                % set the source to the fieldName and the value to the found value.
                results.(cField)(cResults).matches(cResult).source = modelFields{resultEntries(cResult)};
                results.(cField)(cResults).matches(cResult).value = relResults{simorder(cResults),resultEntries(cResult)};
            end
            if printLevel > 0
                % print the individual ids and fields the similarity was achieved on
                fprintf('ID: %s ', model.(cField){relResultPos(simorder(cResults))});
                matchingFields = {results.(cField)(cResults).matches(:).source};
                matchingValues = {results.(cField)(cResults).matches(:).value};
                % filter ID Field
                idpos = strcmp(matchingFields,cField);
                matchingFields = matchingFields(~idpos);
                matchingValues = matchingValues(~idpos);
                if ~isempty(matchingFields)
                   fprintf('with the following matching values:\n');
                   matches = strcat(matchingFields,{': '},matchingValues);
                   fprintf('%s', strjoin(matches,'; '));
                end
                fprintf('\n\n');
            end
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function [LMSout,blms,Rsq]=LMSpol(y,p,x)
%Syntax: [LMSout,blms,Rsq]=LMSpol(y,p,x)
%_______________________________________
%
% Calculates the Least Median of Squares (LMS) polynomial regression
% parameters and output. It searches all the possible combinations
% of points and makes the intercept adjustment for every combination.
%
% LMSout is the LMS estimated values vector.
% blms is the LMS [intercept slopes] vector.
% Rsq is the R-squared.
% y is the vector of the dependent variable.
% p is the order of the polynomial.
% x is the vector of the independent variable.
%
% Reference:
% Rousseeuw PJ, Leroy AM (1987): Robust regression and outlier detection. Wiley.
%
%
% Alexandros Leontitsis
% Institute of Mathematics and Statistics
% University of Kent at Canterbury
% Canterbury
% Kent, CT2 7NF
% U.K.
%
% University e-mail: al10@ukc.ac.uk (until December 2002)
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% Sep 3, 2001.

if nargin<1 | isempty(y)==1
   error('Not enough input arguments.');
else
   % y must be a column vector
   y=y(:);
   % n is the length of the data set
   n=length(y);
end

if nargin<2 | isempty(p)==1
   % If p is omitted give it the value of 1
   p=1;
else
   % p must be a scalar
   if sum(size(p))>2
      error('p must be a scalar.');
   end
   % p must be a non-negative integrer
   if round(p)-p~=0 | p<0
      error('p must be a non-negative integrer');
   end
end

if nargin<3 | isempty(x)==1
   % If x is omitted give it the values 1:n
   x=(1:n)';
else
   % x must be a column vector
   x=x(:);
   % x and y must have the same length
   if n~=size(x,1)
      error('x and y must have the same length.');
   end
end

if n<=p
   error('The ploynomial order is too large for the data set.');
end

% Prepare the matrix X for regression.
for i=1:p
   X(:,i)=x.^i;
end

% Perform the LMS regression
[LMSout, blms, Rsq]=LMSreg(y,X);


text:
function [A] = spm_matrix(P, order)
% returns an affine transformation matrix
% FORMAT [A] = spm_matrix(P, order)
% P(1)  - x translation
% P(2)  - y translation
% P(3)  - z translation
% P(4)  - x rotation about - {pitch} (radians)
% P(5)  - y rotation about - {roll}  (radians)
% P(6)  - z rotation about - {yaw}   (radians)
% P(7)  - x scaling
% P(8)  - y scaling
% P(9)  - z scaling
% P(10) - x affine
% P(11) - y affine
% P(12) - z affine
%
% order (optional) application order of transformations.
%
% A     - affine transformation matrix
%___________________________________________________________________________
%
% spm_matrix returns a matrix defining an orthogonal linear (translation,
% rotation, scaling or affine) transformation given a vector of
% parameters (P).  By default, the transformations are applied in the
% following order (i.e., the opposite to which they are specified):
%
% 1) shear
% 2) scale (zoom)
% 3) rotation - yaw, roll & pitch
% 4) translation
%
% This order can be changed by calling spm_matrix with a string as a
% second argument. This string may contain any valid MATLAB expression
% that returns a 4x4 matrix after evaluation. The special characters 'S',
% 'Z', 'R', 'T' can be used to reference the transformations 1)-4)
% above. The default order is 'T*R*Z*S', as described above.
%
% SPM uses a PRE-multiplication format i.e. Y = A*X where X and Y are 4 x n
% matrices of n coordinates.
%
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_matrix.m 1149 2008-02-14 14:29:04Z volkmar $

% pad P with 'null' parameters
%---------------------------------------------------------------------------
q  = [0 0 0 0 0 0 1 1 1 0 0 0];
P  = [P q((length(P) + 1):12)];

% default multiplication order if not specified
%---------------------------------------------------------------------------
if nargin < 2
    order = 'T*R*Z*S';
end;

T  =   [1   0   0   P(1);
        0   1   0   P(2);
        0   0   1   P(3);
        0   0   0   1];

R1  =  [1    0      0          0;
        0    cos(P(4))  sin(P(4))  0;
        0   -sin(P(4))  cos(P(4))  0;
        0    0      0          1];

R2  =  [cos(P(5))  0    sin(P(5))  0;
        0          1    0      0;
       -sin(P(5))  0    cos(P(5))  0;
        0          0    0          1];

R3  =  [cos(P(6))   sin(P(6))   0  0;
       -sin(P(6))   cos(P(6))   0  0;
        0           0           1  0;
        0           0       0  1];

R   = R1*R2*R3;

Z   =  [P(7)    0       0       0;
        0       P(8)    0       0;
        0       0       P(9)    0;
        0       0       0       1];

S   =  [1       P(10)   P(11)   0;
        0       1   P(12)   0;
        0       0       1   0;
        0       0       0       1];

A = eval(sprintf('%s;', order));
if ~isnumeric(A) || ndims(A) ~= 2 || any(size(A) ~= 4)
    error('Order expression ''%s'' did not return a valid 4x4 matrix.', ...
          order);
end;

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [c, ceq] = findEarliestArrivalNonlconFun(arrivalUT, departUT, dvFunc, maxDV)
    dv = dvFunc(arrivalUT, departUT);
    c(1) = dv - maxDV;
%     c(2) = -(arrivalUT-departUT);

    ceq = [];
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function FbVisualize( FB, show )
% Used to visualize a series of 1D/2D/3D filters.
%
% For 1D and 2D filterabnks also shows the Fourier spectra of the filters.
%
% USAGE
%  FbVisualize( FB, [show] )
%
% INPUTS
%  FB      - filter bank to visualize (either 2D, 3D, or 4D array)
%  show    - [1] figure to use for display
%
% OUTPUTS
%
% EXAMPLE
%  FB=FbMake(1,1,0);  FbVisualize( FB, 1 );  %1D
%  load FbDoG.mat;    FbVisualize( FB, 2 );  %2D
%  FB=FbMake(3,1,0);  FbVisualize( FB, 3 );  %3D
%
% See also FBAPPLY2D, FILTERVISUALIZE
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<2 || isempty(show) ); show=1; end
if( show<=0); return; end;

% get Fourier Spectra for 1D and 2D filterbanks
nd = ndims(FB)-1;
if( nd==1 || nd==2 )
  FBF=zeros(size(FB));
  if( nd==1 )
    for n=1:size(FB,1);  FBF(n,:)=abs(fftshift(fft(FB(n,:)))); end
  else
    for n=1:size(FB,3);  FBF(:,:,n)=abs(fftshift(fft2(FB(:,:,n)))); end
  end
end

% display
figure(show); clf;
if( nd==1 )
  r = (size(FB,2)-1)/2;
  subplot(1,3,1); plot( -r:r, FB );
  subplot(1,3,2); plot( (-r:r)/(2*r+1), FBF );
  subplot(1,3,3); stem( (-r:r)/(2*r+1), max(FBF,[],1) );

elseif( nd==2 )
  subplot(1,3,1); montage2(FB);  title('filter bank');
  subplot(1,3,2); montage2(FBF);
  title('filter bank fft');
  subplot(1,3,3); im(sum(FBF,3));  title('filter bank fft coverage');

elseif( nd==3 )
  n = size(FB,4); nn = ceil( sqrt(n) ); mm = ceil( n/nn );
  for i=1:n
    subplot(nn,mm,i);
    filterVisualize( FB(:,:,:,i), 0 );
  end
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function [csubl, csubd] = aerodyn(alpha)

% aerodynamic properties

% input

%  alpha = angle-of-attack (radians)

% output

%  csubl = lift coefficient (non-dimensional)
%  csubd = drag coefficient (non-dimensional)

% STS aero - J. Betts model

% Orbital Mechanics with Matlab

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rtd = 180.0 / pi;

cl0 = -0.20704d0;

cl1 = 0.029244d0;

cd0 = 0.07854d0;

cd1 = -6.1592d-3;

cd2 = 6.21408d-4;

alphad = rtd * alpha;

csubl = cl0 + cl1 * alphad;

csubd = cd0 + (cd1 + cd2 * alphad) * alphad;


text:
function [d] = vprintf(v,c,varargin)
%
% Copyright Aditya Khosla http://mit.edu/khosla
%
% Please cite this paper if you use this code in your publication:
%   A. Khosla, J. Xiao, A. Torralba, A. Oliva
%   Memorability of Image Regions
%   Advances in Neural Information Processing Systems (NIPS) 2012
%

if(v>=c)
	d=fprintf(varargin{:});
end

author: adikhosla

repo: feature-extraction


text:
% TestScript.m
%
% This script tests the quaternion library functions to ensure that each
% function output is consistent.
%
% Date          Author          Notes
% 27/09/2011    SOH Madgwick    Initial release

%% Start of script

close all;                          % close all figures
clear;                              % clear all variables
clc;                                % clear the command terminal

%% Axis-angle to rotation matrix

axis = [1 2 3];
axis = axis / norm(axis);
angle = pi/2;

R = axisAngle2rotMat(axis, angle);
num = ' % 1.5f';
a = sprintf('\rAxis-angle to rotation matrix:');
b = sprintf(strcat('\r', num, '\t', num, '\t', num), R(1,:));
c = sprintf(strcat('\r', num, '\t', num, '\t', num), R(2,:));
d = sprintf(strcat('\r', num, '\t', num, '\t', num), R(3,:));
disp(strcat(a,b,c,d));

%% Axis-angle to quaternion

q = axisAngle2quatern(axis, angle);
num = ' % 1.5f';
a = sprintf('\rAxis-angle to quaternion:');
b = sprintf(strcat('\r', num, '\t', num, '\t', num, '\t', num), q);
disp(strcat(a,b));

%% Quaternion to rotation matrix

R = quatern2rotMat(q);
num = ' % 1.5f';
a = sprintf('\rQuaternion to rotation matrix:');
b = sprintf(strcat('\r', num, '\t', num, '\t', num), R(1,:));
c = sprintf(strcat('\r', num, '\t', num, '\t', num), R(2,:));
d = sprintf(strcat('\r', num, '\t', num, '\t', num), R(3,:));
disp(strcat(a,b,c,d));

%% Rotation matrix to quaternion

q = rotMat2quatern(R);
num = ' % 1.5f';
a = sprintf('\rRotation matrix to quaternion:');
b = sprintf(strcat('\r', num, '\t', num, '\t', num, '\t', num), q);
disp(strcat(a,b));

%% Rotation matrix to ZYX Euler angles

euler = rotMat2euler(R);
num = ' % 1.5f';
a = sprintf('\rRotation matrix to ZYX Euler angles:');
b = sprintf(strcat('\r', num, '\t', num, '\t', num), euler);
disp(strcat(a,b));

%% Quaternion to ZYX Euler angles

euler = quatern2euler(q);
num = ' % 1.5f';
a = sprintf('\rQuaternion to ZYX Euler angles:');
b = sprintf(strcat('\r', num, '\t', num, '\t', num), euler);
disp(strcat(a,b));

%% ZYX Euler angles to rotation matrix

R = euler2rotMat(euler(1), euler(2), euler(3));
num = ' % 1.5f';
a = sprintf('\rZYX Euler angles to rotation matrix:');
b = sprintf(strcat('\r', num, '\t', num, '\t', num), R(1,:));
c = sprintf(strcat('\r', num, '\t', num, '\t', num), R(2,:));
d = sprintf(strcat('\r', num, '\t', num, '\t', num), R(3,:));
disp(strcat(a,b,c,d));

%% End of script

author: xioTechnologies

repo: Oscillatory-Motion-Tracking-With-x-IMU


text:
function W = randInitializeWeights(L_in, L_out)
    %RANDINITIALIZEWEIGHTS Randomly initialize the weights of a layer with L_in
    %incoming connections and L_out outgoing connections
    %   W = RANDINITIALIZEWEIGHTS(L_in, L_out) randomly initializes the weights
    %   of a layer with L_in incoming connections and L_out outgoing
    %   connections.
    %
    %   Note that W should be set to a matrix of size(L_out, 1 + L_in) as
    %   the column row of W handles the "bias" terms
    %
    EPSILON_INIT = 0.12;
    W = rand(L_out, 1 + L_in) * 2 * EPSILON_INIT - EPSILON_INIT;
end

author: worldveil

repo: coursera-ml


text:
classdef ttopKLogits < matlab.unittest.TestCase
    % ttopKLogits   Unit tests for sampling.topKLogits

    % Copyright 2020 The MathWorks, Inc.

    methods(Test)
        function testForKIsOne(test)
            % When K=1 topKLogits followed by softmax is just the one-hot
            % for the argmax.
            numClasses = 100;
            % Make an arbitrary class the largest logit.
            class = randi([1,numClasses]);
            x = rand([numClasses,1]);
            x(class) = max(x,[],'all')+1;
            % Call topKLogits with K=1
            topK = 1;
            logits = sampling.topKLogits(x,topK);
            logits = dlarray(logits);
            % Apply softmax to get a probability distribution over the
            % classes.
            actProb = softmax(logits,'DataFormat','C');
            % The expected output is a delta distribution on the arbitrary
            % class.
            expProb = zeros([numClasses,1]);
            expProb(class) = 1;
            test.verifyEqual(extractdata(actProb),expProb);
        end

        function testForLargeK(test)
            % When K>1 the topKLogits followed by softmax are only
            % non-zero at the largest K classes. If those classes have
            % equal logits the distribution is uniform.
            numClasses = 100;
            K = 10;
            classes = (1:K)+randi([1,numClasses-K]);
            x = rand([numClasses,1]);
            x(classes) = max(x,[],'all')+1;
            logits = sampling.topKLogits(x,K);
            logits = dlarray(logits);
            actProb = softmax(logits,'DataFormat','C');
            expProb = zeros([numClasses,1]);
            expProb(classes) = 1/K;
            test.verifyEqual(extractdata(actProb),expProb);
        end

        function dlarrayIsSupported(test)
            % Test we can call topKLogits on a dlarray.
            numClasses = 10;
            k = 5;
            x = dlarray(rand([numClasses,1]));
            yAct = sampling.topKLogits(x,k);
            yExp = sampling.topKLogits(extractdata(x),k);
            test.verifyEqual(extractdata(yAct),yExp);
        end
    end

end

author: matlab-deep-learning

repo: transformer-models


text:
% Homogeneous Propagation Medium Example
%
% This example provides a simple demonstration of using k-Wave for the
% simulation and detection of the pressure field generated by an initial
% pressure distribution within a two-dimensional homogeneous propagation
% medium.
%
% author: Bradley Treeby
% date: 29th June 2009
% last update: 17th October 2011
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% k-Wave 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 Lesser General Public License for
% more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

clear all;

% =========================================================================
% SIMULATION
% =========================================================================

% create the computational grid
Nx = 128;           % number of grid points in the x (row) direction
Ny = 128;           % number of grid points in the y (column) direction
dx = 0.1e-3;        % grid point spacing in the x direction [m]
dy = 0.1e-3;        % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);

% define the properties of the propagation medium
medium.sound_speed = 1500;  % [m/s]
medium.alpha_coeff = 0.75;  % [dB/(MHz^y cm)]
medium.alpha_power = 1.5;

% create initial pressure distribution using makeDisc
disc_magnitude = 5; % [Pa]
disc_x_pos = 50;    % [grid points]
disc_y_pos = 50;    % [grid points]
disc_radius = 8;    % [grid points]
disc_1 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

disc_magnitude = 3; % [Pa]
disc_x_pos = 80;    % [grid points]
disc_y_pos = 60;    % [grid points]
disc_radius = 5;    % [grid points]
disc_2 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

source.p0 = disc_1 + disc_2;

% define a centered circular sensor
sensor_radius = 4e-3;   % [m]
num_sensor_points = 50;
sensor.mask = makeCartCircle(sensor_radius, num_sensor_points);

% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor);

% =========================================================================
% VISUALISATION
% =========================================================================

% plot the initial pressure and sensor distribution
figure;
imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, source.p0 + cart2grid(kgrid, sensor.mask), [-1 1]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;

% plot the simulated sensor data
figure;
imagesc(sensor_data, [-1, 1]);
colormap(getColorMap);
ylabel('Sensor Position');
xlabel('Time Step');
colorbar;


text:
function f = rcnn_scale_features(f, feat_norm_mean)
% My initial experiments were conducted on features with an average norm
% very close to 20. Using those features, I determined a good range of SVM
% C values to cross-validate over. Features from different layers end up
% have very different norms. We rescale all features to have an average norm
% of 20 (why 20? simply so that I can use the range of C values found in my
% initial experiments), to make the same search range for C reasonable
% regardless of whether these are pool5, fc6, or fc7 features. This strategy
% seems to work well. In practice, the optimal value for C ends up being the
% same across all features.
target_norm = 20;
f = f .* (target_norm / feat_norm_mean);

author: rbgirshick

repo: rcnn

save_path: github-repos/MATLAB/rbgirshick-rcnn


text:
function [xroot, froot] = brent (f, x1, x2, rtol)

% solve for a single real root of a nonlinear equation

% Brent's method

% input

%  f    = objective function coded as y = f(x)
%  x1   = lower bound of search interval
%  x2   = upper bound of search interval
%  rtol = algorithm convergence criterion

% output

%  xroot  = real root of f(x) = 0
%  froot  = function value at f(x) = 0

% Orbital Mechanics with MATLAB

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global iter;

% machine epsilon

eps = 2.23e-16;

e = 0;

a = x1;
b = x2;

fa = feval(f, a);

fb = feval(f, b);

fc = fb;

for iter = 1:1:50

    if (fb * fc > 0)
       c = a;
       fc = fa;
       d = b - a;
       e = d;
    end

    if (abs(fc) < abs(fb))
       a = b;
       b = c;
       c = a;
       fa = fb;
       fb = fc;
       fc = fa;
    end

    tol1 = 2 * eps * abs(b) + 0.5 * rtol;

    xm = 0.5 * (c - b);

    if (abs(xm) <= tol1 || fb == 0)
       break;
    end

    if (abs(e) >= tol1 && abs(fa) > abs(fb))
       s = fb / fa;

       if (a == c)
          p = 2 * xm * s;
          q = 1 - s;
       else
          q = fa / fc;
          r = fb / fc;
          p = s * (2 * xm * q * (q - r) - (b - a) * (r - 1));
          q = (q - 1) * (r - 1) * (s - 1);
       end

       if (p > 0)
          q = -q;
       end

       p = abs(p);

       min = abs(e * q);

       tmp = 3 * xm * q - abs(tol1 * q);

       if (min < tmp)
          min = tmp;
       end

       if (2 * p < min)
          e = d;
          d = p / q;
       else
          d = xm;
          e = d;
       end
    else
       d = xm;
       e = d;
    end

    a = b;
    fa = fb;

    if (abs(d) > tol1)
       b = b + d;
    else
       b = b + sign(xm) * tol1;
    end

    fb = feval(f, b);

end

xroot = b;
froot = fb;


text:
function mu = moebius ( n )

%*****************************************************************************80
%
%% MOEBIUS returns the value of MU(N), the Moebius function of N.
%
%  Definition:
%
%    MU(N) is defined as follows:
%
%      MU(N) = 1 if N = 1;
%              0 if N is divisible by the square of a prime;
%              (-1)^K, if N is the product of K distinct primes.
%
%  First values:
%
%     N  MU(N)
%
%     1    1
%     2   -1
%     3   -1
%     4    0
%     5   -1
%     6    1
%     7   -1
%     8    0
%     9    0
%    10    1
%    11   -1
%    12    0
%    13   -1
%    14    1
%    15    1
%    16    0
%    17   -1
%    18    0
%    19   -1
%    20    0
%
%    As special cases, MU(N) is -1 if N is a prime, and MU(N) is 0
%    if N is a square, cube, etc.
%
%    The Moebius function is related to Euler's totient function:
%
%      PHI(N) = Sum ( D divides N ) MU(D) * ( N / D ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the value to be analyzed.
%
%    Output, integer MU, the value of MU(N).
%    If N is less than or equal to 0, MU will be returned as -2.
%    If there was not enough internal space for factoring, MU
%    is returned as -3.
%
  if ( n <= 0 )
    mu = -2;
    return
  end

  if ( n == 1 )
    mu = 1;
    return
  end
%
%  Factor N.
%
  [ nfactor, factor, power, nleft ] = i4_factor ( n );

  if ( nleft ~= 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MOEBIUS - Fatal error!\n' );
    fprintf ( 1, '  Not enough factorization space.\n' );
    mu = -3;
    return
  end

  mu = 1;

  for i = 1 : nfactor

    mu = - mu;

    if ( 1 < power(i) )
      mu = 0;
      return
    end

  end

  return
end


text:
function G = greens_function_mono(x,y,z,xs,src,f,conf)
%GREENS_FUNCTION_MONO Green's function in the frequency domain
%
%   Usage: G = greens_function_mono(x,y,z,xs,src,f,conf)
%
%   Input options:
%       x,y,z   - x,y,z points for which the Green's function should be
%                 calculated / m
%       xs      - position of the source
%       src     - source model of the Green's function. Valid models are:
%                   'ps'  - point source
%                   'ls'  - line source
%                   'pw'  - plane wave
%                   'dps' - dipole point source
%       f       - frequency of the source / Hz
%       conf    - configuration struct (see SFS_config)
%
%   Output parameters:
%       G       - Green's function evaluated at the points x,y,z
%
%   GREENS_FUNCTION_MONO(x,y,z,xs,src,f,conf) calculates the Green's function
%   for the given source model located at xs for the given points x,y and the
%   frequency f.
%
%   See also: sound_field_mono

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

%% ===== Checking of input  parameters ==================================
% Disabled checking for performance reasons

%% ===== Configuration ==================================================
c = conf.c;
phase = conf.phase;

%% ===== Computation =====================================================
% Frequency
omega = 2*pi*f;
% Calculate Green's function for the given source model
if strcmp('ps',src)
    % Source model for a point source: 3D Green's function.
    %
    %              1  e^(-i w/c |x-xs|)
    % G(x-xs,w) = --- -----------------
    %             4pi      |x-xs|
    %
    % https://sfs.rtfd.io/en/3.2/sources/#equation-fd-point
    %
    G = 1/(4*pi) * exp(-1i*omega/c .* sqrt((x-xs(1)).^2+(y-xs(2)).^2+(z-xs(3)).^2)) ./ ...
            sqrt((x-xs(1)).^2+(y-xs(2)).^2+(z-xs(3)).^2);

elseif strcmp('dps',src)
    % Source model for a dipole point source: derivative of 3D Green's function.
    %
    %  d                1   / iw       1    \   (x-xs) ns
    % ---- G(x-xs,w) = --- | ----- + ------- | ----------- e^(-i w/c |x-xs|)
    % d ns             4pi  \  c     |x-xs| /   |x-xs|^2
    %
    % r = |x-xs|
    r = sqrt((x-xs(1)).^2+(y-xs(2)).^2+(z-xs(3)).^2);
    % scalar = (x-xs) nxs
    scalar = xs(4).*(x-xs(1)) + xs(5).*(y-xs(2))  + xs(6).*(z-xs(3));
    %
    G = 1/(4*pi) .* (1i*omega/c + 1./r) .* scalar./r.^2 .* exp(-1i*omega/c.*r);

elseif strcmp('ls',src)
    % Source model for a line source: 2D Green's function.
    %
    %                i   (2) / w        \
    % G(x-xs,w) =  - -  H0  |  - |x-xs|  |
    %                4       \ c        /
    %
    % https://sfs.rtfd.io/en/3.2/sources/#equation-fd-line
    %
    G = -1i/4 * besselh(0,2,omega/c* ...
        sqrt( (x-xs(1)).^2 + (y-xs(2)).^2 + (z-xs(3)).^2 ));

elseif strcmp('pw',src)
    % Source model for a plane wave:
    %
    % G(x,w) = e^(-i w/c n x)
    %
    % https://sfs.rtfd.io/en/3.2/sources/#equation-fd-plane
    %
    % Direction of plane wave
    nxs = xs(:,1:3) / norm(xs(:,1:3));
    % Calculate sound field
    G = exp(-1i*omega/c.*(nxs(1).*x+nxs(2).*y+nxs(3).*z));
else
    error('%s: %s is not a valid source model for the Green''s function', ...
        upper(mfilename),src);
end

% Add phase to be able to simulate different time steps
G = G .* exp(-1i*phase);

author: sfstoolbox

repo: sfs-matlab


text:
% Add SDP solver into the MATLAB path

function [success] = initSDP()

%     restoredefaultpath;
     addpath(genpath('../Solver/SDPNAL+v1.0/'), path);
%

end

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function S = oblateSurfaceArea(elli, varargin)
%OBLATESURFACEAREA  Approximated surface area of an oblate ellipsoid.
%
%   S = oblateSurfaceArea(R1,R2)
%
%   Example
%   oblateSurfaceArea
%
%   See also
%   geom3d, ellipsoidSurfaceArea, prolateSurfaceArea
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2015-07-03,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2015 INRA - Cepia Software Platform.

%% Parse input argument

if size(elli, 2) == 7
    R1 = elli(:, 4);
    R2 = elli(:, 5);

elseif size(elli, 2) == 1 && ~isempty(varargin)
    R1 = elli(:, 1);
    R2 = varargin{1};
end

assert(R1 < R2, 'First radius must be smaller than second radius');

% surface theorique d'un ellipsoide oblate
% cf http://fr.wikipedia.org/wiki/Ellipso%C3%AFde_de_r%C3%A9volution
e = sqrt(R2.^2 - R1.^2) ./ R2;
S = 2 * pi * R2.^2 + pi * R1.^2 * log((1 + e) ./ (1 - e)) ./ e;


text:
function a = gk324_inverse ( n, x )

%*****************************************************************************80
%
%% GK324_INVERSE returns the inverse of the GK324 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%   18 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Gregory, David Karney,
%    A Collection of Matrices for Testing Computational Algorithms,
%    Wiley, 1969, page 51,
%    LC: QA263.G68.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real X(N-1), the first N-1 entries of the last row.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n

      if ( i < n )

        if ( j < i )
          a(i,j) = 0.0;
        elseif ( j == i )
          a(i,j) = 1.0 / ( 1.0 - x(i) );
        elseif ( j == i + 1 )
          a(i,j) = - 1.0 / ( 1.0 - x(i) );
        elseif ( i + 1 < j )
          a(i,j) = 0.0;
        end

      elseif ( i == n )

        if ( j == 1 )

          a(i,j) = - x(1) / ( 1.0 - x(1) );

        elseif ( j < n )

          a(i,j) = ( x(j-1) - x(j) ) / ( 1.0 - x(j) ) / ( 1.0 - x(j-1) );

        elseif ( j == n )

          a(i,j) = 1.0 / ( 1.0 - x(n-1) );

        end

      end

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function fx = p15_fun ( option, nvar, x )

%*****************************************************************************80
%
%% P15_FUN evaluates the function for problem 15.
%
%  Title:
%
%    The Trigger Circuit.
%
%  Description:
%
%    The current flow of a trigger circuit with an operational amplifier
%    is modeled.  The variables are voltages, with X(6) the output
%    voltage and X(7) the input voltage.
%
%    The function has the form
%
%      F(X) = A * X + PHI ( X )
%
%    where A is a 6 by 7 matrix, and PHI is a nonlinear term, that is,
%
%      F(I) = SUM ( 1 <= J <= 7 ) A(I,J) * X(J) + PHI ( X )
%
%  Options:
%
%    Melhem lists the following limit points in X(7):
%
%    ( 0.04936  0.54735  0.04944  0.04944  0.12920  1.16602  0.60185 )
%    ( 0.23577  0.66296  0.23759  0.23760  0.62083  9.60913  0.32286 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 August 20008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Rami Melhem, Werner Rheinboldt,
%    A Comparison of Methods for Determining Turning Points of Nonlinear Equations,
%    Computing,
%    Volume 29, Number 3, September 1982, pages 201-226.
%
%    Gerd Poenisch, Hubert Schwetlick,
%    Computing Turning Points of Curves Implicitly Defined by Nonlinear
%    Equations Depending on a Parameter,
%    Computing,
%    Volume 26, Number 2, June 1981, pages 107-121.
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Input, integer NVAR, the number of variables.
%
%    Input, real X(NVAR), the argument of the function.
%
%    Output, real FX(NVAR-1), the value of the function at X.
%

%
%  Get the linear coefficients.
%
  array = p15_gx ( );
%
%  Compute the linear portion of the function.
%
  fx(1:nvar-1) = 0.0;

  for i = 1 : nvar - 1
    for j = 1 : nvar
      fx(i) = fx(i) + array(i,j) * x(j);
    end
  end
%
%  Add the nonlinear terms.
%
  fx(2) = fx(2) + 5.6D-08 * ( exp ( 25.0 * x(2) ) - 1.0 );
  fx(5) = fx(5) + 5.6D-08 * ( exp ( 25.0 * x(5) ) - 1.0 );
  fx(6) = fx(6) - 7.65 * atan ( 1962.0 * ( x(3) - x(1) ) ) / 0.201;

  return
end


text:
classdef MimRemapCyan < MimGuiPlugin
    % MimRemapCyan. Gui Plugin for setting paint colour
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimRemapCyan is a Gui Plugin for the MIM Toolkit.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Right lower'
        SelectedText = 'Right lower'
        ToolTip = 'Marks airway as right lower'
        Category = 'Airway label'
        Visibility = 'Dataset'
        Mode = 'Edit'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 5
        ButtonHeight = 1
        Icon = 'paint.png'
        IconColour = GemMarkerPoint.DefaultColours{4}
        Location = 34
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.PaintBrushColour = 4;
        end

        function enabled = IsEnabled(gui_app)
             enabled = gui_app.IsDatasetLoaded && gui_app.ImagePanel.OverlayImage.ImageExists && ...
                isequal(gui_app.ImagePanel.SelectedControl, 'Map');
        end

        function is_selected = IsSelected(gui_app)
            is_selected = gui_app.ImagePanel.PaintBrushColour == 4;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function test_old_trialdef

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY

% This script tests the implementation of the new representation of trial
% specific information in the data structure.
% The idea is as follows:
% -instead of relying on data.cfg.previous....trl, the trial information
% (when consistent with data representation, so only on the level of raw
% data, and in original sampling rate) should be directly on the main level
% of the structure.
% -this means that raw data should get a field sampleinfo containing the
% first 2 columns of the original trl-matrix.
% -the 3rd column of the trl matrix is the offset column and is represented
% in each trial's individual time axis.
% -user-defined columns > 3 should go to the field trialinfo. this field
% can percolate deeper into the pipeline. yet, once data is not raw
% datatype anymore (or resampled) the sampleinfo field should be removed,
% because consistency with the data is lost.
% -the sampleinfo and trialinfo fields will be updated by subselecting rows
% -sampleinfo should be changed by ft_redefinetrial
% -sampleinfo should be removed by ft_resampledata and any other function
% working on the raw datatype and outputting any other datatype
% -sampleinfo and trialinfo should be concatenated in ft_appenddata
% -sampleinfo and trialinfo should be adjusted by ft_rejectartifact

% needed for the dccnpath function, since we will change directory later on
addpath(fileparts(mfilename('fullpath')));

cd(dccnpath('/home/common/matlab/fieldtrip/data/test/original/meg/ctf275/A0132_Aud-Obj-Recognition_20051115_02.ds'));
headerfile = 'A0132_Aud-Obj-Recognition_20051115_02.res4';
datafile   = 'A0132_Aud-Obj-Recognition_20051115_02.meg4';
%hdr        = ft_read_header(headerfile);

cfg          = [];
cfg.datafile = datafile;
cfg.trl      = [[1001:1000:10001]' [2000:1000:11000]' round(randn(10,1)*100)];
cfg.trl(:,4) = [ones(5,1); ones(5,1)*2];
cfg.continuous = 'yes';
data1          = ft_preprocessing(cfg);
cfg.trl(:,1:2) = cfg.trl(:,1:2) + 11000;
cfg.trl(:,3)   = round(randn(10,1)*100);
data2          = ft_preprocessing(cfg);

% test subselection of trialsclear
cfg        = [];
cfg.trials = (1:5);
datax      = ft_preprocessing(cfg, data1);

% test some other preprocessing (so nargin==2 for ft_preprocessing)
cfg        = [];
cfg.demean = 'yes';
datax2     = ft_preprocessing(cfg, data1);

% test ft_checkdata
datay = data1;
datay = rmfield(datay,   'sampleinfo');
datay = rmfield(datay,   'trialinfo');
datay = ft_checkdata(datay, 'hastrialdef', 'yes');

datay.trialinfo(:) = 3;
datay = ft_checkdata(datay, 'hastrialdef', 'yes'); % should give all(datay.trialinfo==3)

datay.trialinfo = [];
datay = ft_checkdata(datay, 'hastrialdef', 'yes'); % should give all(datay.trialinfo==trl(:,4))

datay = data1;
datay = rmfield(datay,   'sampleinfo');
datay = rmfield(datay,   'trialinfo');
datay.cfg = rmfield(datay.cfg, 'trl');
datay = ft_checkdata(datay, 'hastrialdef', 'yes');

% test appenddata
data1b = rmfield(data1, 'sampleinfo');
data1b = rmfield(data1b, 'trialinfo');
%data1b.cfg = rmfield(data1b.cfg, 'trl');
data2b = rmfield(data2, 'sampleinfo');
data2b = rmfield(data2b, 'trialinfo');
%data2b.cfg = rmfield(data2b.cfg, 'trl');
data3 = ft_appenddata([], data1b, data2b);

% test ft_redefinetrial with and without sampleinfo and trialinfo present in data
cfg = [];
cfg.toilim = [-0.5 0.5];
data1bx    = ft_redefinetrial(cfg, data1b);
data1x     = ft_redefinetrial(cfg, data1 );

cfg = [];
cfg.begsample =    1 + round(rand(10,1)*200);
cfg.endsample = 1000 - round(rand(10,1)*200);
data1by    = ft_redefinetrial(cfg, data1b);
data1y     = ft_redefinetrial(cfg, data1 );

cfg = [];
cfg.trl      = data1.cfg.trl;
cfg.trl(:,1) = cfg.trl(:,1)+round(rand(10,1)*20);
cfg.trl(:,2) = cfg.trl(:,2)-round(rand(10,1)*20);
data1bz      = ft_redefinetrial(cfg, data1b);
data1z       = ft_redefinetrial(cfg, data1 );

% test ft_resampledata
cfg = [];
cfg.resamplefs = 300;
cfg.detrend    = 'no';
cfg.demean     = 'yes';
data1rs        = ft_resampledata(cfg, data1);
data1brs       = ft_resampledata(cfg, data1b);
cfg.trials     = [1:5];
data1rs2       = ft_resampledata(cfg, data1b);

% test ft_databrowser
cfg = [];
cfg.viewmethod      = 'butterfly';
cfg.preproc.demean  = 'yes';
cfg.channel         = 'MEG';
cfg.continuous      = 'yes';
ft_databrowser(cfg, data1b);

% test ft_rejectartifact
cd(dccnpath('/home/common/matlab/fieldtrip/data/test/original/meg/ctf275/A0132_Aud-Obj-Recognition_20051115_02.ds'));
headerfile = 'A0132_Aud-Obj-Recognition_20051115_02.res4';
datafile   = 'A0132_Aud-Obj-Recognition_20051115_02.meg4';
hdr        = ft_read_header(headerfile);

cfg          = [];
cfg.datafile = datafile;
cfg.trl      = [[1001:1000:10001]' [2000:1000:11000]' round(randn(10,1)*100)];
cfg.trl(:,4) = [ones(5,1); ones(5,1)*2];
cfg.continuous = 'yes';
cfg.artfctdef.type  = 'eog';
cfg.artfctdef.eog.channel = 'MLO11';
cfg.artfctdef.eog.feedback = 'no'; % yes does not work in the test script, see http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=2840
cfg          = ft_rejectartifact(cfg);

datay = data1;
datay = rmfield(datay,   'sampleinfo');
datay = rmfield(datay,   'trialinfo');
datay2       = ft_rejectartifact(cfg, datay);

% test all other functions using raw data as input (to remove sampleinfo)

% test all high level functions to correctly pass on the trialinfo field
% ft_timelockanalysis
cfg              = [];
cfg.vartrllength = 2;
timelock         = ft_timelockanalysis(cfg, data1);

%ft_denoise_synthetic
cfg          = [];
cfg.gradient = 'G1BR';
dataG1BR     = ft_denoise_synthetic(cfg, data1);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function maxcorrel=plotcorrelatioondistance(delta_vec,dir1,dir2, name1, name2, marker, col, savethis)
% maxcorrel=plotcorrelatioondistance(delta_vec,res_dir, res_name, marker,
% color savethis)

method='average';
ii=1;
for delta = delta_vec

    for iteration =1:10;
        res_dir = [dir1 num2str(delta*100) dir2];
        res_name = [name1 num2str(iteration) name2];
        load ([res_dir '/' res_name]);
%         load ([peval.data_path peval.data_dir '/' peval.data_name]); %data
        [Wxk,Hkt,centers,Vxkpix]=reshapeGaP(res.hvec,res.cxcy,peval);
        Vxtpixbg=reshape(Wxk*Hkt,peval.nx,peval.ny,peval.nt)+peval.bg;
        resid=(Vxtpixbg-res.dpixc);
        resid_norm=resid./sqrt(Vxtpixbg);
        %         [Z,H,T,perm] = dendrogram_subtreepixels(resid_norm,'average', p , centers(:,1)+1, centers(:,2)+1, savethis)
        data=resid_norm;
        sized = size(data);
        dveccr= reshape(data,sized(1)*sized(2), sized(3));
        ccd = (corrcoef(dveccr'));
%         ccds=squareform(1-ccd);
        e=eye(size(ccd));
        ccdflip=e-(ccd-e);
        ccds=squareform(1-ccdflip);

        Z = linkage(ccds,method);

        if savethis
            save ([res_dir '/corrcoefdist.mat'], 'Z')
        end
        maxcorrel(ii,iteration)=1-min(Z(:,3));
    end
    ii=ii+1;
end

% style='o--r';
% marker='o';
% col='r';
style = [marker col '--'];
figure(1)
hold on
deltanm=delta_vec*106;
m=mean(maxcorrel,2);
s=std(maxcorrel,[],2);
mincorr=min(maxcorrel,[],2);
errorbar(deltanm, m,s,style)
plot(deltanm,mincorr,['-v' col], 'linewidth',2)
xlabel('Separation of sources [nm]')
ylabel('Maximum correlation in residuals')
grid on

figure(2)
hold on
m=mean(maxcorrel,2);
s=std(maxcorrel,[],2);
errorbar(delta_vec, m,s,style)
xlabel('Separation of sources [pix]')
plot(delta_vec,mincorr,['-v' col], 'linewidth',2)
ylabel('Maximum correlation in residuals')
grid on

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function class=fixTopology(classFile,datFile,outFile)
%
% class=fixTopology(classFile,datFile,outFile)
%
% Function: Given a mrGray .Class file and the corresponding unsegmented
% data from a .dat file, fixes the topology and writes the result to a new
% .Class file.
%
%
% Requires:
% removeBridgesMex (in Matlab path)
%
% HISTORY:
%   110202 ISS wrote it.
%   2002.11.13 RFD (bob@white.stanford.edu) cleaned up and updated to use
%   the new bridgeRemovalMex file.
%
%

warndlg('fixTopology has been replaced by mrgRemoveBridges.  Calling that program for you.')
mrgRemoveBridges(classFile,datFile,outFile);
return;

if ~exist('classFile','var')
    [f,p]=uigetfile('*.class;*.Class', 'Select a class file...');
    if(isnumeric(f) & f==0) return; end
    classFile = fullfile(p,f);
    disp(['Class file: ',classFile]);
end
if ~exist('datFile','var')
    [f,p]=uigetfile('*.dat;*.DAT', 'Select corresponding vAnatomy...');
    if(isnumeric(f) & f==0) return; end
    datFile = fullfile(p,f);
    disp(['Anatomy file: ',datFile]);
end
if(~exist('outFile','var') | isempty(outFile))
    [f,p]=uiputfile('*.class;*.Class', 'When finished, save Class file as...');
    if(isnumeric(f) & f==0) return; end
    outFile = fullfile(p,f);
    disp(['Output Class file: ',outFile]);
end

%-----------------------Read in presegmented data----------------------------
% Data is coming in permuted!
preSegClass = readClassFile(classFile,0,0);
preSegWhite = preSegClass.data;
preSegWhite(preSegClass.data~=16)=0;
% The handle removal function sets all non-zero values to white matter, so
% we can leave the 16's in there. We just need to set all non-white matter
% to zero.

%-----------------------Read in unsegmented data-----------------------------

[imageData,mmPerVox,img_dim]=readVolAnat(datFile);
% [ysize,xsize,zsize]=size(imageData);
% imageData=reshape(imageData,[ysize,xsize,zsize]);

% Permuting to match preseg data
imageData=permute(imageData,[2 1 3]);
scale = 225/256;
imageData = uint8(double(imageData)*scale);

%-----------------------Call removeBridgesMex----------------------------

if(isfield(preSegClass.header,'params') & length(preSegClass.header.params)>=6)
    % the 6 classification params are csf mean, gray mean, white mean,
    % noise stdev, confidence and smoothness.
    avgWhite = preSegClass.header.params(3)*scale;
    avgGray = preSegClass.header.params(2)*scale;
    threshold = mean([avgWhite,avgGray]);
else
    warning('Classification parameters not found in class file header- making some guesses...');
    avgWhite = mean(imageData(preSegWhite==16));
    threshold = 0.9*avgWhite;
    avgGray = 0.8*avgWhite;
end
tic;
corrected = removeBridgesMex(preSegWhite,imageData,avgWhite,threshold,avgGray);
toc;
% Voxels marked with 236 should be removed (make a cut).
% Voxels marked with 245 should be added.
% Voxels marked with 226 are 'alternative' voxels (what are they?).
% Voxels marked with 225 are unchanged voxels.
%preSegClass.data(corrected==236)=0;
%preSegClass.data(corrected==245 | corrected==225) = 16;

% For some reason, this works better than the previous, even though they
% *should* be equivalent. I suspect that the bridge removal code doesn't
% properly mark all the removed voxels. However, we need to put the CSF
% back in.
preSegWhite = uint8(corrected==225 | corrected==245);
preSegWhite(preSegWhite>0) = preSegClass.type.white;
preSegWhite(preSegClass.data==preSegClass.type.csf) = preSegClass.type.csf;
preSegClass.data = preSegWhite;

class = writeClassFile(preSegClass, outFile);

% class = writeClassFileFromRaw(preSegClass.data, outFile, [0,0; 16,16; 32,32; 48,48]);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ y2, m2, d2, f2 ] = ymdf_prev_roman ( y1, m1, d1, f1 )

%*****************************************************************************80
%
%% YMDF_PREV_ROMAN returns the Roman YMDF date of the previous day.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1, M1, D1, real F1,
%    the YMDF date.
%
%    Output, integer Y2, M2, D2, real F2,
%    yesterday's YMDF date.
%
  y2 = y1;
  m2 = m1;
  d2 = d1 - 1;
  f2 = f1;

  [ y2, m2, d2 ] = day_borrow_roman ( y2, m2, d2 );

  return
end


text:
function [out] = prep_grandAverage(dat)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prep_grandAverage
%
% Synopsis:
%   [out] = prep_grandAverage(dat,<var>)
%
% Example :
%   [out] = prep_grandAverage({dat1,dat2,...,datN});
%
% Arguments:
%     dat - Segmented data, structure or data itself, in a cell type
%
% Returns:
%     out - Data structure averaged across subjects
%
% Description:
%    This function averages the data across subjects. It is recommended that
%    input data be separated by each class in advance.
%     continuous data should be [time * channels]
%     epoched data should be [time * channels * trials]
%
% See also 'https://github.com/PatternRecognition/OpenBMI'
%
% Min-ho Lee, 01-2018
% mh_lee@korea.ac.kr
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s = length(dat);
x = cell(1,s);
if isstruct(dat{1})
    if ~prod(cellfun(@(x) isfield(x,'x'),dat))
        warning('OpenBMI: Data structure must have a field named ''x''');return
    end
    EachDataSize=cellfun(@(y) size(y.x), dat,'UniformOutput',false);
    if ~isequal(EachDataSize{:})
        warning('OpenBMI: Data size must be same (Epoched or continuous)');return
    end
% chan, class °°ÀºÁö È®ÀÎ ÇÊ¿ä.
    for i=1:s
        x{i} = dat{i}.x;
    end
    nd = ndims(dat{1}.x);
elseif isnumeric(dat{1}) && (ismatrix(dat{1}) || ndims(dat{1})==3)
    EachDataSize=cellfun(@size, dat,'UniformOutput',false);
    if ~isequal(EachDataSize{:})
        warning('OpenBMI: Data size must be same (Epoched or continuous)');return
    end
    x = dat;
    nd = ndims(dat{1});
else
    warning('OpenBMI: Check for format or dimension of the data');return
end

switch nd
    case 2
        Dat = cat(3,x{:});
        M = mean(Dat,3);
    case 3
        Dat = cat(4,x{:});
        M = mean(Dat,4);
end

if isstruct(dat{1})
    out = rmfield(dat{1},'x');
    out.x = M;
elseif isnumeric(dat{1})
    out = M;
end

author: PatternRecognition

repo: OpenBMI


text:
function M_sh2pv = beamWeightsPressureVelocity(basisType)
% BEAMWEIGHTSPRESSUREVELOCITY Convert from SH signals to pressure velocity
%
%   Returns the 4x4 matrix that converts first-order SH signals to pressure
%   velocity signals, with the ordering [p v_x v_y v_z], where v_xyz are the
%   cartesian components of the acoustic veclocity vector.
%
%   basisType:  'complex' or 'real' for the respective type of spherical
%       harmonics under use.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% BEAMWEIGHTSPRESSUREVELOCITY.M - 11/7/2013
% Archontis Politis, archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch basisType
    case 'real'
        M_sh2pv = sqrt(4*pi)*[  1   0           0           0;
                                0   0           0           1/sqrt(3);
                                0   1/sqrt(3)   0           0;
                                0   0           1/sqrt(3)   0];
    case 'complex'
        M_sh2pv = sqrt(4*pi)*[  1   0           0           0;
                                0   1/sqrt(6)   0           -1/sqrt(6);
                                0   -1i/sqrt(6) 0           -1i/sqrt(6);
                                0   0           1/sqrt(3)   0];
end

author: polarch

repo: Spherical-Array-Processing


text:
function X = rocholBackSub(ch, Y);

% ROCHOLBACKSUB Backsubstitute the representation of the rank one Cholesky.

% ROCHOL

%/~
% This would be the long way of doing it.
%L = rocholExtract(ch);
%X1 = L'\Y;
%~/
X = zeros(size(Y));
t = zeros(1, size(Y, 2));
X(ch.n, :) = Y(ch.n, :)/ch.s(ch.n);
for i = ch.n-1:-1:1
  t = t + ch.v(i+1)*X(i+1, :);
  X(i, :) = Y(i, :)/ch.s(i) - ch.u(i)*t;
end

%/~
%disp(max(max(X - X1)));
%~/

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat