text:
%% imx
% Below is a demonstration of the features of the |imx| function

%%
clear; close all; clc;

%% Syntax
% |hf=imx(varargin);|

%% Description
% The imx (image explorer) function is a GUI to navigate and segment 3D
% image data.

%% Examples
%
%% Example: Segmenting MRI data

%%
% Get a 3D image (see also |dcmFolder2MATobject| to import DICOM data)
testCase=2;

defaultFolder = fileparts(fileparts(mfilename('fullpath'))); %Set main folder
loadSavePath=fullfile(defaultFolder,'data','imseg'); %Path for loading and saving segmentation data
switch testCase
    case 1 %MATLAB brain data
        load mri;
        M=squeeze(D); %example image data set
        v=2./[1,1,.4]; %example voxel size
    case 2 %MRI imported from DICOM files, see also HELP_dcmFolder2MATobject
        pathName=fullfile(defaultFolder,'data','DICOM','0001_human_calf');
        loadName=fullfile(pathName,'IMDAT','IMDAT.mat');

        IMDAT_struct=load(loadName); %The image data structure
        G = IMDAT_struct.G; %Geometric/spatial information
        v = G.v; %The voxel size
        M = IMDAT_struct.type_1; %The image data
end

%%
% Start segmentation using |imx|
hf=imx(M,v,loadSavePath);

%%
% *Viewing*
%
% * _Toolbar_ - Note all functions displayed in the graphical user
% interface toolbar, equipt with an information panel and descriptions
% of keyboard shortcuts just below the toolbar.
% Instructions for each widget will appear along this top bar when in use.
% Right click to exit from any given widget.
%
% * _View manipulation_ - Use the View Control Widget to zoom, pan, and
% rotate a figure using different mouse gestures (right mouse button
% for zoom, left for panning, and middle/scroll for rotating).
% Use the sliding bars on the left hand side of the graphical user
% interface to change slice view from all 3 axes.
% Press home to return to active slice.
% Use the colorbar button to change the limits of the color axis displayed
% below the slice, shifting the color contrast across a slice and
% increasing the visibility of contour lines.
%
% <<HELP_imx_View.gif>>
%
% *Sketching a contour*
%
% * _Sample_ - Begin segmenting an element with s = sample sketch contour.
% With the activated widget, click on an area of contrast that outlines the
% element to generate an approximate dotted contour. If the generated outline is
% unsuitable, click elsewhere along the contrast to produce a new contour.
% Once the sample contour is acceptable, right click out of the
% widget. See below for an example of this process in segmenting the skin.
%
% <<HELP_imx_Sample.gif>>
%
% * _Cut_ - Use c = cut sketched contour to edit a sample sketch by
% boxing and removing an incorrect section of the outline.
%
% <<HELP_imx_Cut.gif>>
%
% * _Draw_ - After cutting a portion of a sample sketch, use the drawing
% tool to complete the sketch if needed. Use d = draw contour to
% place connected sketch points along a contour where the path needs to be filled in.
% After exiting the widget, observe that the drawn sketch will be in a
% different color compared to the original.
%
% <<HELP_imx_Draw.gif>>
%
% * _Adjust_ - To edit the alignment of certain points along a bad contour,
% use the adjust tool. To move a sketch point to a new location, click on
% the point and then click on its new placement.
%
% <<HELP_imx_Adjust.gif>>
%
% * _Delete_ - Use the delete tool to erase any sketch contour on the
% slice.
%
% <<HELP_imx_Delete.gif>>
%
%
% *Accepting a contour*
%
% * _Select_ - Once a sketch contour is well fit to an element on the given
% slice, use a = accept contour to select the contour, changing it to a
% solid green path. If multiple sketches exist on the slice, select your
% desired option from the displayed window to either add a new contour,
% merge with another contour, or replace the existing accepted contours.
% See below for an example of merging multiple contours of interest.
%
% <<HELP_imx_Accept.gif>>
%
% * _Convert_ - Use the convert button to easily switch from an accepted
% contour back to a sketch contour if further adjustments are needed. Use
% the reset button to return to the original, unsegmented slice.
%
% <<HELP_imx_Convert.gif>>
%
% * _Grow/Shrink_ - Use + = grow contour and - = shrink contour to
% manipulate the size of an accepted contour to best fit a given element.
% Growing and shrinking contour lines will likely wrinkle the path and
% require smoothing (see Smooth).
%
% <<HELP_imx_GrowShrink.gif>>
%
% * _Smooth_ - Use q = smooth accepted contour to smoothen a contour path.
% If more extreme smoothing is necessary, press p while the smoothing tool
% is activated to adjust the smoothing parameter accordingly.
%
% <<HELP_imx_Smooth.gif>>
%
% *Finalizing segmentation* - After the slice is complete, press space to
% go to the next slice and continue segmentation for all MRI slices.
% Periodically check the view of all slices using the showHide tool to
% verify continuity along the length of an element until completion.
%
% <<HELP_imx_ShowHide.gif>>
%

%%
%
% <<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:
% Buttons for loading/saving a prgram from commandLine
function interface_save_load_CommandLine_program_panel(h,xPos,yPos,loadSave_commandLine_Panel_width,...
    loadSave_commandLine_Panel_height)
%% About:
% Construct the interface-panel for load/save a program from commandLine
% interface

%% Areguments:
% h: handle of the main frame
% frame_width: width of the main frame (pixels)
% frame_height: height of the main frame (pixels)

% Copyright: Mohammad SAFEEA, 20th-July-2018

panel_Width=loadSave_commandLine_Panel_width;
panel_height=loadSave_commandLine_Panel_height;

x0=xPos;
y0=yPos;

pnl_h = uipanel(h,'Title','My Programs',...
                'Units','pixels');
            %     'FontSize',12,...
            %    'BackgroundColor','white',...
pnl_h.Position=[x0 y0 panel_Width panel_height];

% number of buttons
I=2;
J=1;
% separation ratio
rh=0.1;
rv=0.3;
btns_width=panel_Width/(J+(J+1)*rh);
btns_height=panel_height/(I+(I+1)*rv);

texts={'Load Program';
        'Save program'};

callback_functions={'fcnCallback_cmd_Load_coammnd_line_Program';
    'fcnCallback_cmd_save_coammnd_line_Program'};

icon_image={'icon_load.jpeg';
    'icon_save.jpeg'};
handles_Array=[];
for i=1:I
    for j=1:J
        labels_text=texts{i,j};
        btn_handle= uicontrol(pnl_h,'Style', 'pushbutton', 'String', labels_text,...
                'Units','pixels',...
                'Callback', callback_functions{i,j});
        x0=panel_Width-(btns_width*(1+rh))*(J+1-j);
        y0=panel_height-(btns_height*(1+rv))*(i)-0.5*rv*btns_height;
        btn_handle.HorizontalAlignment='center';
        btn_handle.Position=[x0 y0 btns_width btns_height];
        btn_handle.FontWeight='bold';
        % Add an icon to the button
        cd icons;
        d1=0;
        if btns_width<btns_height
            d1=btns_width;
        else
            d1=btns_height;
        end
        imageName=icon_image{i,j};
        addImageToButton(btn_handle,imageName,d1);
    end
end

end

function addImageToButton(btn_handle,imageName,d1)
        im=imread(imageName);
        scaleFactor=d1/max(max(max(size(im))));
        width=scaleFactor*size(im,2);
        width=floor(width);
        height=scaleFactor*size(im,1);
        height=floor(height);
        im2 = imresize(im,[width,height]);
        set(btn_handle,'CData',im2);
        cd ..;
end

repo: KST-Kuka-Sunrise-Toolbox


author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
% f6.m
% Schaffer's F6 function
% commonly used to test optimization/global minimization problems
%
% z = 0.5+ (sin^2(sqrt(x^2+y^2))-0.5)/((1+0.01*(x^2+y^2))^2)

function [out]=f6(in)
 x=in(:,1);
 y=in(:,2);
 num=sin(sqrt(x.^2+y.^2)).^2 - 0.5;
 den=(1.0+0.01*(x.^2+y.^2)).^2;

 out=0.5 +num./den;

author: HuangCongQing


text:
function[h,hl]=wavespecplot(varargin)
%WAVESPECPLOT  Plot of wavelet spectra together with time series.
%
%   WAVESPECPLOT(T,X,P,W) where T is the time, X is a time series, and W is
%   the wavelet transform at periods P, makes a two-component plot.  The
%   upper subplot has the time series X plotted against its time axis T.
%   The lower subplot has the transform W, or its modulus ABS(W) if W is
%   complex-valued, plotted versus the time axis T and the periods P.
%
%   WAVESPECPLOT(T,X,P,W,R) optionally plots ABS(W).^R instead of W.
%
%   WAVESPECPLOT(T,X,P,W,R,CI) makes a filled contour plot with contour
%   intervals CI.  If CI is not input, then the spectrum is plotted using
%   PCOLOR, which is faster to render but slow to print. For making final
%   figures, it is better to use the contouring option.
%
%   WAVESPECPLOT(T,X,P,W1,W2,...WN,...) makes an N+1 component plot, with
%   W1 in the second subplot, W2 in the third, etc.
%
%   After plotting, all subplots are packed together using PACKFIG.
%   H=WAVESPECPLOT(...) returns the handles to the subplots.
%
%   [H,HL]=WAVESPECPLOT(...) also returns the handles to the lines plotted
%   in the first panel.
%
%   If X is complex-valued, both the real and imaginary parts are plotted
%   in the uppermost subplot.
%
%   Usage: h=wavespecplot(t,x,p,w);
%          h=wavespecplot(t,x,p,w,r);
%          [h,hl]=wavespecplot(t,x,p,w,r,ci);
%
%   'wavespecplot --f' generates a sample figure.
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2004--2015 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1},'--f')
   type makefigs_wavespecplot;
   makefigs_wavespecplot
   return
end

t=varargin{1};
u=varargin{2};
p=varargin{3};
t=t(:);
p=p(:);

nspec=nargin-3;

n=1;
ci=[];

if length(varargin{end})==1
  if length(varargin{end-1})==1
     ci=varargin{end};
     n=varargin{end-1};
     nspec=nspec-2;
  else
     n=varargin{end};
     nspec=nspec-1;
  end
else
  if isvector(varargin{end})
     ci=varargin{end};
     n=varargin{end-1};
     nspec=nspec-2;
  end
end

a=min(t);
b=max(t);

if isreal(u)
  c=maxmax(abs(u));
else
  c=max([maxmax(abs(real(u))) maxmax(abs(imag(u)))]);
end
c=c*1.1;

subplot(nspec+1,1,1)
if ~isreal(u)
  [hl1,hl2]=uvplot(t,u);
  hl=[hl1;hl2];
else
  hl=plot(t,u);
end
xlim([a,b])
ylim([-c c])
grid off
%hlines(0,'k:')

for i=1:nspec
  subplot(nspec+1,1,i+1)
  if ~isreal(varargin{3+i}) || n~=1
        varargin{3+i}=abs(varargin{3+i});
  end
  T=(varargin{3+i}').^n;
  if isempty(ci)
    pcolor(t,p,T),shading interp
  else
    contourf(oprod(1+0*p,t),oprod(p,1+0*t),T,ci),hold on
    caxis([min(ci) max(ci)])
  end
  xlim([a,b])
  if p(1)<p(end)
      flipy
  end
  hold on
  set(gca,'tickdir','out')
  ylog
end
h=packfig(nspec+1,1,'rows');

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
 function subplot_stack(x, ys, str_title, colors)
%function subplot_stack(x, ys, str_title, colors)
% a tight stack of subplots to show L signal components
% in
%	x	[N,1]
%	ys	[N,L]		or ?

if nargin == 1 && streq(x, 'test'), subplot_stack_test, return, end
if nargin < 2, ir_usage, end
if ~isvar('colors') || isempty(colors), colors = {'c', 'y'}; end
if ~isvar('str_title') || isempty(str_title)
	str_title = '';
end

if ~iscell(ys)
	ys = {ys};
end
L = size(ys{1},2);

apos = get(gca, 'position'); % current axes position
for ll=1:L
%	pos = [0.1 0.1+0.8/L*(L-ll) 0.8 0.8/L];
	pos = [apos(1) apos(2)+apos(4)/L*(L-ll) apos(3) apos(4)/L];
	subplot('position', pos) % l b w h
	for ip=1:length(ys)
		plot(	x, real(ys{ip}(:,ll)), colors{1+2*(ip-1)}, ...
			x, imag(ys{ip}(:,ll)), colors{2+2*(ip-1)})
		if ip == 1, hold on, end
	end
	hold off
	axis tight
	ytick(0), set(gca, 'yticklabel', '')
	fontsize = 10;
	fontweight = 'normal';
	texts(1.02, 0.8, sprintf('%d.', ll), ...
		'fontsize', fontsize, 'fontweight', fontweight)

	if ll==1, title(str_title), end
	if ll<L
		xtick off
	end
end

function subplot_stack_test
x = linspace(0,1,101)';
y = exp(2i*pi*x*[1:5]);
if im
	clf, subplot(121)
	subplot_stack(x, y)
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function conv2pool()
    global config mem;
    curr_layer_idx = config.misc.current_layer;
    mem.layer_inputs{curr_layer_idx} = reshape(accumarray(mem.pooling_matrix{curr_layer_idx}, mem.activations{curr_layer_idx-1}(:), ...
                        [size(mem.activations{curr_layer_idx-1},1)*size(mem.activations{curr_layer_idx-1},2)/4, 1]), size(mem.activations{curr_layer_idx-1}, 2)/4, size(mem.activations{curr_layer_idx-1}, 1));
end

author: jimmy-ren


text:
% STD_CHANGROUP - Create channel groups for plotting.
%
% Usage:
%                >> STUDY = std_changroup(STUDY, ALLEEG);
%                >> STUDY = std_changroup(STUDY, ALLEEG, chanlocs, 'interp');
% Inputs:
%   ALLEEG     - Top-level EEGLAB vector of loaded EEG structures for the dataset(s)
%                in the STUDY. ALLEEG for a STUDY set is typically loaded using
%                POP_LOADSTUDY, or in creating a new STUDY, using POP_CREATESTUDY.
%   STUDY      - EEGLAB STUDY set comprising some or all of the EEG datasets in ALLEEG.
%   chanlocs   - EEGLAB channel structure. Only construct the STUDY.changrp
%                structure for a subset of channels.
%   'interp'   - optional input in case channel locations are interpolated
%
% Outputs:
%   STUDY      - The input STUDY set structure modified according to specified user
%                edits, if any. The STUDY.changrp structure is created. It contains as
%                many elements as there are channels. For example, STUDY.changrp(1)
%                is the first channel. Fields of the changrp structure created at this
%                point are
%                    STUDY.changrp.name      : name of the channel group
%                    STUDY.changrp.channels  : cell array containing channel labels
%                                              for the group.
%                    STUDY.changrp.setinds   : indices of datasets containing the
%                                              selected channels.
%                    STUDY.changrp.allinds   : indices of channels within the datasets
%                                              above.
%
% Authors: Arnaud Delorme, CERCO, 2006

% Copyright (C) Arnaud Delorme, CERCO, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function STUDY = std_changroup(STUDY, ALLEEG, alllocs, interp);

if nargin < 4
    interp = 'off';
end

% union of all channel structures
% -------------------------------
inputloc = 0;
if nargin >= 3
    if ~isempty(alllocs)
        inputloc = 1;
    end
end
if ~inputloc
    alllocs = eeg_mergelocs(ALLEEG.chanlocs);
end

% create group for each electrode
% -------------------------------
if isstruct(alllocs)
    alllocs = { alllocs.labels };
end
STUDY.changrp = [];
for indc = 1:length(alllocs)
    STUDY.changrp(indc).name     = alllocs{indc};
    STUDY.changrp(indc).channels = { alllocs{indc} };
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
% Derive  -- Equations of Motion - 2D (planar) quad-rotor
%
% Model:
%   We will assume that the quadrotor is modeled by two point-masses that
%   are connected by a rigid mass-less bar. There are two motors, each of
%   which apply a force orthogonal to the connecting bar. Assume that air
%   drag is negligable.
%

clc; clear;

%%%% Create symbolic variables:
syms x y q 'real'  % horizontal and vertical position, angle (+i axis = 0)
syms dx dy dq 'real'  % rates
syms ddx ddy ddq 'real' % accelerations
syms u1 u2 'real'  % force applied by rotors
syms d m g 'real'  % distance from center to rotor, rotor mass, gravity

%%%% Unit vectors:

% Inertial reference frame:
i = sym([1;0;0]);
j = sym([0;1;0]);
k = sym([0;0;1]);

% Body reference frame:
e = cos(q)*i + sin(q)*j;  % Unit vector from center to rotor 2
n = -sin(q)*i + cos(q)*j;  % Unit vector orthogonal to e

%%%% Position vectors:
Pc = x*i + y*j;   % Position of the center of the quad rotor
P1 = Pc - d*e;  % Position of rotor 1
P2 = Pc + d*e;  % Position of rotor 2

%%%% Define the state vector and its derivative
z = [x;y;q;dx;dy;dq];  %State
dz = [dx;dy;dq;ddx;ddy;ddq];  %Derivative of state

%%%% Kinematics:

% Define the derivative operator. CHAIN RULE.
derivative = @(f)( jacobian(f,z)*dz );

% Velocities
dPc = derivative(Pc);
dP1 = derivative(P1);
dP2 = derivative(P2);

% Accelerations
ddPc = derivative(dPc);
ddP1 = derivative(dP1);
ddP2 = derivative(dP2);

%%%% Define all force vectors:
Fg1 = -m*g*j;  %Force of gravity on mass 1
Fg2 = -m*g*j;  %Force of gravity on mass 2
Fu1 = u1*n;  %Force due to rotor 1
Fu2 = u2*n;  %Force due to rotor 1

%%%% Force Balance:
sumForces = Fg1 + Fg2 + Fu1 + Fu2;
massAccel = m*ddPc;
eqn1 = dot(sumForces-massAccel,i);
eqn2 = dot(sumForces-massAccel,j);

%%%% Angular Momentum Balance:   (about center of quad rotor)
sumTorques = ...
    cross(P1-Pc, Fg1) + ...
    cross(P2-Pc, Fg2) + ...
    cross(P1-Pc, Fu1) + ...
    cross(P2-Pc, Fu2);
angMomentum = ...
    cross(P1-Pc, m*ddP1) + ...
    cross(P2-Pc, m*ddP2);
eqn3 = simplify(dot(sumTorques-angMomentum, k));

%%%% Collect and solve equations:
vars = [ddx;ddy;ddq];  %This is what we want to find  (the accelerations)
eqns = [eqn1;eqn2;eqn3];  %These are the dynamics equations
[M,f] = equationsToMatrix(eqns, vars); % EoM are linear in acceleration

% In this problem, the equations of motion are simple, so we can solve them
% analytically using the "\" command. For more complicated systems, it is
% typically much faster to do this solve numerically at run-time.
soln = simplify(M\f);
ddxSoln = soln(1);
ddySoln = soln(2);
ddqSoln = soln(3);

%%%% Automatically write out dynamics function:
% This will automatically optimize the function to minimize computation
% time. Not too important here, but makes a big difference for complicated
% functions. It also prevents copy-paste errors.
matlabFunction(...
    ddxSoln, ddySoln, ddqSoln,...
    'file','autoGen_dynamics.m',...
    'vars',{q,u1,u2,m,g,d},...
    'outputs',{'ddx','ddy','ddq'});

%%%% Display the solution to the user:
disp(['ddx = ' ddxSoln]);
disp(['ddy = ' ddySoln]);
disp(['ddq = ' ddqSoln]);

author: MatthewPeterKelly

repo: OptimTraj


text:
function[]=makefigs_morlfreq
%MAKEFIGS_MORLFREQ  Makes a sample figure for MORLFREQ.

fm=linspace(0,4,1000);
[fn,fmin]=morlfreq(fm);

figure
plot([fn' fn'],[fm' fmin']),axis equal,
dlines(1,'k--')
axis([0 4 -1 4]),%vlines(1,'k:'),
hlines(0,'k:')

xlabel('Carrier Frequency (Radian)')
ylabel('Frequency of Max and Min (Radian)')
title('Frequencies of Morlet Maximum and Minimum')


text:
function [doseArray, isCompress, isRemote] = getDoseArray(doseIndex, planC)
%"getDoseArray"
%   Returns the doseArray stored in a planC, in slot doseIndex.  The last
%   returned doseArray is cached in a persistent variable to avoid delays
%   in loading remotely stored or compressed doses.
%
%   If a dose array is compressed it is decompressed before being returned.
%
%   If planC is not passed in, the global planC is used.
%
%   If doseIndex is not valid, or if no valid planC can be found an error is
%   returned.
%
%   isCompress and isRemote are booleans to let the calling function know
%   if the doseArray was compressed and/or remote.
%
%   JRA 10/14/04
%
%Usage:
%   function doseArray = getDoseArray(doseIndex, planC)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

persistent lastDoseStruct;
persistent lastDoseArray;
persistent isLastDoseCompressed;
persistent isLastDoseRemote;

%Check if a doseStruct, not a doseIndex was passed.
if isstruct(doseIndex) && isfield(doseIndex, 'doseArray')

    doseStruct = doseIndex(1);

else
    %An index was passed, extract the doseStruct.
	if ~exist('planC', 'var') %wy
        global planC
	end

	if ~iscell(planC)
        error('Cannot get dose. PlanC is not a valid cell array.');
	end

	indexS = planC{end};

	if doseIndex > length(planC{indexS.dose}) | doseIndex < 1
        error('Cannot get dose.  Requested doseIndex does not exist.')
	end

    doseStruct = planC{indexS.dose}(doseIndex);

end

doseArray = doseStruct.doseArray;

%If remote or compressed, and struct is same as last time, return cached array.
if (isCompressed(doseArray) || ~isLocal(doseArray)) && isequal(doseStruct, lastDoseStruct);
    doseArray = lastDoseArray;
    isCompress = isLastDoseCompressed;
    isRemote = isLastDoseRemote;
    return;
%If remote or compressed and NOT the same as last time, clear cache.
elseif (isCompressed(doseArray) || ~isLocal(doseArray)) && ~isequal(doseStruct, lastDoseStruct);
	lastDoseArray           = [];
	lastDoseStruct          = [];
	isLastDoseCompressed    = [];
	isLastDoseRemote        = [];
end

%To be output for calling function's use.
isCompress  = 0;
isRemote    = 0;

%Decompress and follow all file pointers until get to an array.
while isCompressed(doseArray) || ~isLocal(doseArray);
    if isCompressed(doseArray)
        doseArray   = decompress(doseArray);
        isCompress  = 1;
    elseif ~isLocal(doseArray);
        doseArray   = getRemoteVariable(doseArray);
        isRemote    = 1;
    end
end

lastDoseArray           = doseArray;
lastDoseStruct          = doseStruct;
isLastDoseCompressed    = isCompress;
isLastDoseRemote        = isRemote;

%temp

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function T=findBaylissSubarrays(xyPoints,numLevels,sidelobedB,N,adjMat)
%%FINDBAYLISSSUBARRAYS Given a elements for a linear array or for a planar
%          array with a circular aperture, group the elements into
%          subarrays that discretize Bayliss difference patterns to a given
%          number of levels. For a planar array, both horizontal and
%          vertical difference pattern are overlapped so that the subarrays
%          can produce difference patterns in either direction. This
%          implements the basic method of determining subarrays as in [1].
%
%INPUTS: xyPoints A 2XnumPoints set of numPoints points in the aperture
%             plane corresponding to element positions. For a linear array,
%             this is a 1XnumPoints set of points. The points should be
%             shifted such that the center of the aperture is the origin.
%             The width of the aperture is the maximum distance of any
%             point from the center.
%   numLevels The number of levels to use in discretizing the Bayliss
%             pattern. If this parameter is omitted or an empty matrix is
%             passed, a default of 4 is used.
%  sidelobedB The number of decibels of the ratio of the close-in
%             sidelobe voltages to the main lobe voltage in the differnce
%             pattern. This must be a negative number. A typical value is
%             -30.
%           N The Bayliss tapering is computed using a certain number of
%             terms. Using too many terms can be undesirable as edge
%             illumination increases, as noted in [1]. If this parameter is
%             omitted or an empty matrix is passed, then the default of 17
%             is used. In [1], it is suggested that N be chosen to be
%             <2*a/lambda, where a is the radius of the aperture and
%             lambda the wavelength.
%      adjMat This parameter is only used with planar arrays and is
%             optional. This is a numPointsXnumPoints adjacency matrix such
%             that adjMat(i,j) is nonzero if element i is adjacent to
%             element j. The value of adjMat(i,i) does not matter. Only the
%             lower half of this matrix is used. If this matrix is omitted
%             or an empty matrix is passed, all elements are taken to be
%             adjacent to each other. The purpose of this matrix is to make
%             sure that no disjoint subarrays are formed. That is, a
%             subarray consists of a continuum of adjacent elements with no
%             breaks.
%
%OUTPUTS: T A numSubarraysXnumPoints boolean matrix (a matrix such that
%           T(i,:) is a set of boolean values indicating which elements are
%           in each subarray). The subarrays do not overlap.
%
%The design of subarrays is a very difficult optimization problem. This
%implements the algorithm of [1], which is based on the idea that the
%ability to form good difference beams should play a role in subarray
%design. The method presented is significatnly simpler than many of the
%genetic optimization methods in the literature.
%
%EXAMPLE 1:
%This is an example of a linear array being broken into six subarrays.
% N=17;
% sidelobedB=-30;
% Nx=41;%There are 2*Nx+1 points total.
% %Generate points symmetric about the origin.
% xPoints=(-(Nx-1)/2):1/2:((Nx-1)/2);
% numLevels=6;%Discretize into four levels.
% T=findBaylissSubarrays(xPoints,numLevels,sidelobedB,N);
%
% %Display the elements in the array with different coloring for each
% %subarray.
% figure(1)
% clf
% hold on
% els=xPoints(T(1,:));
% numEls=length(els);
% scatter(els,zeros(1,numEls),'or')
% els=xPoints(T(2,:));
% numEls=length(els);
% scatter(els,zeros(1,numEls),'xg')
% els=xPoints(T(3,:));
% numEls=length(els);
% scatter(els,zeros(1,numEls),'+b')
% els=xPoints(T(4,:));
% numEls=length(els);
% scatter(els,zeros(1,numEls),'*k')
% els=xPoints(T(5,:));
% numEls=length(els);
% scatter(els,zeros(1,numEls),'sc')
% els=xPoints(T(6,:));
% numEls=length(els);
% scatter(els,zeros(1,numEls),'dm')
%
%EXAMPLE 2:
%This is an example of a circular planar array. This requires the formation
%of two Bayliss patterns that are them discretized and overlapped. The
%elements in one row are offset by lambda/4 from those in the next, so they
%do not form a very regular grid.
% %First, we create a circular array. The element locations are given in
% %terms of the wavelength lambda, so lambda will not appear in the
% %equations for the sum beam.
% [xyVals,adjMat]=getShaped2DLattice([49;49],'circular');
%
% sidelobedB=-30;
% N=5;
% numLevels=3;
% T=findBaylissSubarrays(xyVals,numLevels,sidelobedB,N,adjMat);
%
% %Display the subarrays with different colors and symbols.
% figure(2)
% clf
% hold on
% axis square
% numSubarrays=size(T,1);
% dispOpts={'.b','og','xr','+c','*m','sy','dk','vg','^b','<r','>c','pm','hy'};
% for curSubarray=1:numSubarrays
%     points=xyVals(:,T(curSubarray,:));
%     optVal=mod(curSubarray,length(dispOpts))+1;
%     scatter(points(1,:),points(2,:),dispOpts{optVal},'linewidth',2);
% end
% h1=xlabel('x');
% h2=ylabel('y');
% title('Coloring Represents Subarrays')
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%REFERENCES:
%[1] U. R. O. Nickel, "Subarray configurations for digital beamforming with
%    low sidelobes and adaptive interference suppression," in Record of the
%    IEEE International Radar Conference, Alexandria, VA, 8-11 May 1995,
%    pp. 714-719.
%
%August 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(N))
    N=17;
end

if(nargin<3||isempty(sidelobedB))
    sidelobedB=-30;
end

if(nargin<2||isempty(numLevels))
   numLevels=4;
end

switch(size(xyPoints,1))
    case 1%Linear array
        %Bayliss weights are all imaginary.
        g=BaylissLinearTapering(sidelobedB,N,xyPoints);
        T=discretizeAndOverlapVals(imag(g),[],numLevels);
    case 2%Planar array
        numEls=size(xyPoints,2);

        if(nargin<5||isempty(adjMat))
            %If no adjacency matrix is given, then everything is adjacent.
            adjMat=ones(numEls,numEls);
        end

        %Bayliss weights are all complex.
        %Weights for horizontal differencing (in u)
        gHoriz=BaylissTapering(sidelobedB,N,xyPoints);
        %Weights for vertical differencing (in v)
        gVert=BaylissTapering(sidelobedB,N,flipud(xyPoints));

        %Bayliss weights are all imaginary.
        T=discretizeAndOverlapVals(imag(gHoriz),imag(gVert),numLevels,adjMat);
    otherwise
        error('Only linear and planar arrays are implemented.')
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function write_tsurf_to_gltf(filename, tsh)
  drawnow;
  % Check that plot is set up for per-vertex data
  assert(strcmp(tsh.FaceColor,'interp'));
  % Even though you may pass it  with 'CData' seems it gets moved to
  % 'FaceVertexCData'
  assert(size(tsh.FaceVertexCData,2) == 1);
  assert(size(tsh.FaceVertexCData,1) == size(tsh.Vertices,1));

  % get color map
  CM = tsh.Parent.Colormap;
  CA = caxis(tsh.Parent);
  % pad colormap to power of 2
  pad = 2^ceil(log2(size(CM,1))) - size(CM,1);
  CA = [CA(1) CA(1) + (CA(2)-CA(1)) * (size(CM,1)+pad)/size(CM,1)];
  CM(end+(1:pad),:) = repmat(CM(end,:),pad,1);
  im = permute(flipud(CM),[1 3 2]);

  D = (tsh.FaceVertexCData-CA(1))/(CA(2)-CA(1));
  R = axisangle2matrix([1 0 0],pi/2);
  VN = -tsh.VertexNormals;
  if ~isempty(VN);
   VN = VN*R;
  end

  writeGLTF_and_validate( ...
    filename, ...
    tsh.Vertices * R , ...
    tsh.Faces, ...
    'Normals',VN, ...
    'MagFilter','NEAREST', ...
    'MinFilter','NEAREST', ...
    'TextureCoordinates',[D D], ...
    'TextureImage',im);
end

author: alecjacobson

repo: gptoolbox


text:
function [grad]= B_linear(future_layers, curr_layer)
grad = GetFutureGrad(future_layers, curr_layer);
end

author: singaxiong

repo: SignalGraph


text:
function linplus_test62 ( )

%*****************************************************************************80
%
%% TEST62 tests R8VM_DET.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST62\n' );
  fprintf ( 1, '  R8VM_DET, determinant of a Vandermonde matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
%
%  Set the matrix.
%
  [ a, seed ] = r8vm_random ( n, n, seed );

  r8vm_print ( n, n, a, '  The Vandermonde matrix:' );
%
%  Copy the matrix into a general array.
%
  a2 = r8vm_to_r8ge ( n, n, a );
%
%  Compute the determinant.
%
  det = r8vm_det ( n, a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  R8VM_DET computes the determinant = %14.6e\n', det );
%
%  Factor the general matrix.
%
  [ a2_lu, pivot, info ] = r8ge_fa ( n, a2 );
%
%  Compute the determinant.
%
  det = r8ge_det ( n, a2_lu, pivot );

  fprintf ( 1, '  R8GE_DET computes the determinant = %14.6e\n', det );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function test_corenums

end

author: ckczzj

repo: CHAN

save_path: github-repos/MATLAB/ckczzj-CHAN


text:
function haveNewerRelease = tapas_compare_versions(online_version,offline_version)
    %% Check if online_version is higher.
    %
    % Input
    %       online_version      -- Online version as string of numbers separated
    %                               by dots, i.e. '3.3.0'
    %       offline_version     -- Online version as string of numbers separated
    %                               by dots, i.e. '2.7.0.3'
    %
    % Output
    %       haveNewerRelease    -- true if version number of online_version is
    %                               higher than that of offline_version

    % muellmat@ethz.ch
    % copyright (C) 2020
    %

    haveNewerRelease = false; % If both releases are the same, return false
    online_numbers = str2double(strsplit(online_version,'.'));
    offline_numbers = str2double(strsplit(offline_version,'.'));
    % Check whether they have the same lenght (oterwise fill with zeros)
    n_online = numel(online_numbers);
    n_offline = numel(offline_numbers);
    if n_online > n_offline
        offline_numbers(n_online) = 0; % matlab fills rest with zeros
    elseif n_online < n_offline
        online_numbers(n_offline) = 0; % matlab fills rest with zeros
    end
    % Now that they have the same length, compare them number by number.
    %   Return at the first difference.
    for ind = 1:max(n_online,n_offline)
        if online_numbers(ind) > offline_numbers(ind)
            haveNewerRelease = true;
            return;
        elseif online_numbers(ind) < offline_numbers(ind)
            haveNewerRelease = false;
            return;
        end
    end

end

author: translationalneuromodeling

repo: tapas


text:
function [bcx,bcy] = specific_flow(xbd,ybd)
%tightcavity_flow   Reference problem 5.3 alternative inflow condition
%   [bcx,bcy] = specific_flow(xbd,ybd);
%   input
%          xbd          x coordinate vector
%          ybd          y coordinate vector
%
%   specifies watertight cavity flow boundary condition
%   IFISS function: DJS; 6 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
bcx=0*xbd; bcy=0*xbd;
k=find(ybd==1 & xbd>-1 & xbd<1); bcx(k)=1;
return


text:
function Torques(nlp, sys)

    domains = sys.Gamma.Nodes.Domain;

    for i=1:numel(domains)
        domain = domains{i};

        u = domain.Inputs.Control.u;
        u2r = tovector(norm(u).^2);
        u2r_fun = SymFunction(['torque_' sys.Gamma.Nodes.Name{i}],u2r,{u});
        rs_phase = getPhaseIndex(nlp,sys.Gamma.Nodes.Name{i});
        addRunningCost(nlp.Phase(rs_phase),u2r_fun,{'u'});

    end
    nlp.update;
end

author: ayonga

repo: frost-dev

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


text:
function test_suite=test_assert_true()
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_assert_true_exceptions

    assertExceptionThrown(@()assertTrue(...
                                false),...
                                        'assertTrue:falseCondition');
    assertExceptionThrown(@()assertTrue(...
                                struct),...
                                        'assertTrue:invalidCondition');
    assertExceptionThrown(@()assertTrue(...
                                [true true]),...
                                        'assertTrue:invalidCondition');

function test_assert_true_passes
    assertTrue(true);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [fLogLikelihood] = calc_AkiLikelihood(mCatalog, fBValue, fBinning)
% function [fLogLikelihood] = calc_AkiLikelihood(mCatalog, fBValue, fBinning)
% ---------------------------------------------------------------------------
% Calculates the likelihood of a b-value fit
%
% Input parameters:
%   mCatalog        Earthquake catalog
%   fBValue         b-value
%   fBinning        Binning of the earthquake magnitudes (default 0.1)
%
% Output parameters:
%   fLogLikelihood  Log-likelihood
%
%@ARTICLE{Aki1965,
%  author =       "K. Aki",
%  title =        "Maximum likelihood estimate of $b$ in the formula
%                  $\log N = a-bM$ and its confidence limits",
%  journal =      "Bull. Earthquake Re. Inst., Tokyo Univ.",
%  year =         "1965",
%  volume =       "43",
%  pages =        "237-239",
%}
%
% Copyright (C) by Danijel Schorlemmer
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the
% Free Software Foundation, Inc.,
% 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

if ~exist('fBinning', 'var')
  fBinning = 0.1
end

fBPrime = fBValue/(log10(exp(1)));
fMinMag = min(mCatalog(:,6))-(fBinning/2);

fL = ones(length(mCatalog(:,1)),1)*nan;

fL = log(fBPrime) - (fBPrime * (mCatalog(:,6) - fMinMag));
fLogLikelihood = sum(fL);

author: CelsoReyes


text:
% TEXT Add text descriptions to data points
%    TEXT(x,y,str) adds a text description to one or more data points in the
%    current axes using the text specified by str. To add text to one point,
%    specify x and y as scalars in data units. To add text to multiple
%    points, specify x and y as vectors with equal length.
%
%    TEXT(x,y,z,str) positions the text in 3-D coordinates.
%
%    TEXT(...,Name,Value) specifies text properties using one or more
%    Name,Value pair arguments. For example, 'FontSize',14 sets the font
%    size to 14 points. You can specify text properties with any of the
%    input argument combinations in the previous syntaxes. If you specify
%    the Position and String properties as Name,Value pairs, then you do not
%    need to specify the x, y, z, and str inputs.
%
%    TEXT(container,...) creates the text in the axes, group, or transform
%    specified by container, instead of in the current axes.
%
%    T = TEXT(...) returns one or more text objects. Use T to modify
%    properties of the text objects after they are created. For a list of
%    properties and descriptions, see Text Properties. You can specify an
%    output with any of the previous syntaxes.
%
%    Execute GET(T), where T is a text object, to see a list of text object
%    properties and their current values.
%    Execute SET(T) to see a list of text object properties and legal
%    property values.
%
%    See also XLABEL, YLABEL, ZLABEL, TITLE, GTEXT, LINE, PATCH.
%
%    Reference page in Doc Center
%       doc text
%
%

author: jmaih


text:
% Test file for sing/compose.m

function pass = test_compose(pref)

if ( nargin < 1 )
    pref = chebfunpref();
end

% Generate a few random points to use as test values.
seedRNG(6178);
d = 2;
x = 2*(1-10^(-d)) * rand(100, 1) - (1-10^(-d));

% Composition of two SINGFUNs - OP(F, G):
data.exponents = [-1 0];
f = singfun(@(x) sin(x)./(x + 1), data, pref);
g = singfun(@(x) cos(x)./(x + 1), data, pref);
data.exponents = [-1 0];
h = compose(f, @plus, g, data, pref);
op = @(x) (sin(x) + cos(x))./(x + 1);
hVals = feval(h, x);
hExact = op(x);
err = hVals - hExact;
pass(1) = norm(err, inf) < 1e3*eps*get(h, 'vscale');

% Composition of an operator and a SINGFUN - OP(F)
data.exponents = [0.5 0];
f = singfun(@(x) sqrt(x + 1), data, pref);
h = compose(f, @sin);
op = @(x) sin(sqrt(x+1));
hVals = feval(h, x);
hExact = op(x);
err = hVals - hExact;
pass(2) = norm(err, inf) < 1e1*eps*get(h, 'vscale');

% Composition of a SMOOTHFUN and a SINGFUN - G(F)
data.exponents = [0.5 0];
f = singfun(@(x) sqrt(x + 1), data, pref);
g = chebtech2(@(x) cos(x));
h = compose(f, g);
op = @(x) cos(sqrt(x + 1));
hVals = feval(h, x);
hExact = op(x);
err = hVals - hExact;
pass(3) = norm(err, inf) < 1e4*eps*get(h, 'vscale');

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function ExampleMusicSpeechClassification(cDatasetPath)

    if (nargin<1)
        % this script is written for the GTZAN music/speech dataset
        % modify this path or use the function parameter to specify your
        % dataset path
        cDatasetPath = 'd:\dataset\music_speech\';
    end
    if (exist('ComputeFeature') ~= 2)
        error('Please add the ACA scripts (https://github.com/alexanderlerch/ACA-Code) to your path!');
    end
    if ((exist([cDatasetPath 'music']) ~= 7) || (exist([cDatasetPath 'speech']) ~= 7))
        error('Dataset path wrong or does not contain music/speech folders!')
    end

    iNumFeatures = 2;

    % read directory contents
    music_files = dir([cDatasetPath 'music/*.au']);
    speech_files = dir([cDatasetPath 'speech/*.au']);

    v_music = zeros(iNumFeatures,size(music_files,1));
    v_speech = zeros(iNumFeatures,size(speech_files,1));

    % extract features, this may take a while...
    for i = 1:size(music_files, 1)
        v_music(:, i) = ExtractFeaturesFromFile_I(...
            [cDatasetPath 'music/' music_files(i).name]);
    end
    for i = 1:size(speech_files, 1)
        v_speech(:, i) = ExtractFeaturesFromFile_I(...
            [cDatasetPath 'speech/' speech_files(i).name]);
    end

    % assign class labels for training and eval
    C = [zeros(1, size(music_files, 1)) ones(1, size(speech_files, 1))];

    % normalize features
    v = [v_music, v_speech];
    m = mean(v, 2);
    s = std(v, 0, 2);
    v = (v - repmat(m, 1, size(music_files, 1) + size(speech_files, 1)))./...
        repmat(s, 1, size(music_files, 1)+size(speech_files, 1));

    % compute the overall accuracy with cross validation
    [acc, mat] = ToolLooCrossVal(v, C);

    disp('confusion matrix:'),
    disp(mat);

    disp('micro accuracy:'),
    disp(sum(diag(mat)) / sum(sum(mat)))
    tmp = zeros(size(mat, 1), 1);
    for i = 1:size(mat, 1)
        tmp(i) = mat(i, i) / sum(mat(i, :));
    end
    disp('macro accuracy:'),
    disp(mean(tmp))

    % compute the individual feature performance
    [acc1, mat1] = ToolLooCrossVal(v(1, :), C);
    sprintf('centroid accuracy: %f', acc1)
    [acc2, mat2] = ToolLooCrossVal(v(2 ,:), C);
    sprintf('rms accuracy: %f', acc2)
end

function [v] = ExtractFeaturesFromFile_I(cFilePath)

    cFeatureNames = char('SpectralCentroid',...
    'TimeRms');

    % read audio
    [x, fs] = audioread(cFilePath);
    x = x / max(abs(x));

    % compute first feature
    feature = ComputeFeature (deblank(cFeatureNames(1, :)), x, fs);
    v(1, 1)  = mean(feature);

    % compute second feature
    feature = ComputeFeature (deblank(cFeatureNames(2, :)), x, fs);
    v(2, 1)  = std(feature(1,:));
end

author: alexanderlerch

repo: ACA-Code


text:
clear all;
n=400000;
p=1000;
density=0.01;

% generate random data
format compact;
randn('seed',0);
rand('seed',0);
X=sprandn(p,n,density);
nrm=sqrt(sum(X.^2));
X=X/mean(nrm);
%mean_nrm=mean();
%X=X/mean_nrm;

% generate some true model
z=double(sign(full(sprandn(p,1,0.05))));
y=X'*z;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EXPERIMENT 1: Lasso
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('EXPERIMENT FOR LASSO\n');
nrm=sqrt(sum(y.^2));
y=y+0.01*nrm*randn(n,1);    % add noise to the model
nrm=sqrt(sum(y.^2));
y=y*(sqrt(n)/nrm);

clear param;
param.regul='l1';        % many other regularization functions are available
param.loss='square';     % only square and log are available
param.numThreads=-1;    % uses all possible cores
param.normalized=false;  % if the columns of X have unit norm, set to true.
param.averaging_mode=0;  % no averaging, averaging was not really useful in experiments
param.weighting_mode=0;  % weights are in O(1/sqrt(n))  (see help mexStochasticProx)
param.optimized_solver=true;  % best is not to touch this option
param.verbose=false;

% set grid of lambda
max_lambda=max(abs(X*y))/n;
tablambda=max_lambda*(2^(1/4)).^(0:-1:-20);  % order from large to small
param.lambda=tablambda;    % best to order from large to small
tabepochs=[1 2 3 5 10];  % in this script, we compare the results obtained when varying the number of passes over the data.

%% The problem which will be solved is
%%   min_beta  1/(2n) ||y-X' beta||_2^2 + lambda ||beta||_1
fprintf('EXPERIMENT: ALL LAMBDAS IN PARALLEL\n');
% we try different experiments when varying the number of epochs.
% the problems for different lambdas are solve INDEPENDENTLY in parallel
obj=[];
objav=[];
for ii=1:length(tabepochs)
   param.iters=tabepochs(ii)*n;
   fprintf('EXP WITH %d PASS\n',tabepochs(ii));
   nlambdas=length(param.lambda);
   Beta0=zeros(p,nlambdas);
   tic
   [Beta tmp]=mexStochasticProx(y,X,Beta0,param);
   toc
   yR=repmat(y,[1 nlambdas]);
   fprintf('Objective functions: \n');
   obj=[obj; 0.5*sum((yR-X'*Beta).^2)/n+param.lambda.*sum(abs(Beta))];
   obj
   if param.averaging_mode
      objav=[objav; 0.5*sum((yR-X'*tmp).^2)/n+param.lambda.*sum(abs(tmp))];
      objav
   end
   fprintf('Sparsity: \n');
   sum(Beta ~= 0)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EXPERIMENT 2: L2 logistic regression
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('EXPERIMENT FOR LOGISTIC REGRESSION + l2\n');
y=sign(y);
param.regul='l2';
param.loss='logistic';
obj=[];
objav=[];
for ii=1:length(tabepochs)
   param.iters=tabepochs(ii)*n;   % one pass over the data
   fprintf('EXP WITH %d PASS\n',tabepochs(ii));
   nlambdas=length(param.lambda);
   Beta0=zeros(p,nlambdas);
   tic
   [Beta tmp]=mexStochasticProx(y,X,Beta0,param);
   toc
   yR=repmat(y,[1 nlambdas]);
   fprintf('Objective functions: \n');
   obj=[obj;sum(log(1.0+exp(-yR .* (X'*Beta))))/n+0.5*param.lambda.*sum(abs(Beta.^2))];
   obj
   if param.averaging_mode
      objav=[objav; sum(log(1.0+exp(-yR .* (X'*tmp))))/n+0.5*param.lambda.*sum(abs(tmp))];
      objav
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EXPERIMENT 3: L1 logistic regression
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('EXPERIMENT FOR LOGISTIC REGRESSION + l1\n');
y=sign(y);
param.regul='l1';        % many other regularization functions are available
param.loss='logistic';     % only square and log are available
obj=[];
objav=[];
for ii=1:length(tabepochs)
   param.iters=tabepochs(ii)*n;   % one pass over the data
   fprintf('EXP WITH %d PASS\n',tabepochs(ii));
   nlambdas=length(param.lambda);
   Beta0=zeros(p,nlambdas);
   tic
   [Beta tmp]=mexStochasticProx(y,X,Beta0,param);
   toc
   yR=repmat(y,[1 nlambdas]);
   fprintf('Objective functions: \n');
   obj=[obj; sum(log(1.0+exp(-yR .* (X'*Beta))))/n+param.lambda.*sum(abs(Beta))];
   obj
   if param.averaging_mode
      objav=[objav; sum(log(1.0+exp(-yR .* (X'*tmp))))/n+param.lambda.*sum(abs(tmp))];
      objav
   end
   fprintf('Sparsity: \n');
   sum(Beta ~= 0)
end


text:
function varargout = spm_api_nmm(varargin)
% SPM_API_NMM M-file for spm_api_nmm.fig
%      SPM_API_NMM, by itself, creates a new SPM_API_NMM or raises the existing
%      singleton*.
%
%      H = SPM_API_NMM returns the handle to a new SPM_API_NMM or the handle to
%      the existing singleton*.
%
%      SPM_API_NMM('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in SPM_API_NMM.M with the given input arguments.
%
%      SPM_API_NMM('Property','Value',...) creates a new SPM_API_NMM or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before spm_api_nmm_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to spm_api_nmm_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help spm_api_nmm

% Last Modified by GUIDE v2.5 17-Oct-2008 16:14:27

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
    'gui_Singleton',  gui_Singleton, ...
    'gui_OpeningFcn', @spm_api_nmm_OpeningFcn, ...
    'gui_OutputFcn',  @spm_api_nmm_OutputFcn, ...
    'gui_LayoutFcn',  [] , ...
    'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before spm_api_nmm is made visible.
%--------------------------------------------------------------------------
function spm_api_nmm_OpeningFcn(hObject, eventdata, handles, varargin)
set(gcf,'name','Neural mass modelling: review of priors');
try
    handles.DCM = varargin{1};
    unpack_Callback(hObject, eventdata, handles);
end

% Choose default command line output for spm_api_nmm
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes spm_api_nmm wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = spm_api_nmm_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

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

% --- Executes on button press in load.
%--------------------------------------------------------------------------
function load_Callback(hObject, eventdata, handles)

[f,p]       = uigetfile('*.mat','please select DCM file'); cd(p)
name        = fullfile(p,f);
DCM         = load(name,'-mat');
DCM         = DCM.DCM;
DCM.name    = name;
handles.DCM = DCM;
set(handles.name,'String',f);

% display priors
%--------------------------------------------------------------------------
unpack_Callback(hObject, eventdata, handles);
guidata(hObject,handles);

% --- Executes on button press in save.
%--------------------------------------------------------------------------
function save_Callback(hObject, eventdata, handles)
try
    [p,file] = fileparts(handles.DCM.name);
catch
    try
        [p,file] = fileparts(handles.DCM.xY.Dfile);
        file     = ['DCM_' file];
    catch
        file     = ['DCM'  date];
    end
end
[file,fpath]     = uiputfile(['DCM*.mat'],'DCM file to save',file);

if fpath
    handles.DCM.name = fullfile(fpath,file);
    set(handles.name,'String',file);
    DCM              = handles.DCM;
    save(DCM.name,'DCM', spm_get_defaults('mat.format'))
    cd(fpath)
end

% assign in base
%--------------------------------------------------------------------------
assignin('base','DCM',handles.DCM)
guidata(hObject,handles);

% --- Executes on button press in kernels.
%--------------------------------------------------------------------------
function kernels_Callback(hObject, eventdata, handles)
clear spm_erp_L spm_gen_erp
DCM   = handles.DCM;
pst   = str2num(get(handles.pst,'String'));
Hz    = str2num(get(handles.Hz, 'String'));

% initialise states
%--------------------------------------------------------------------------
M     = DCM.M;
model = DCM.options.model;
[x f] = spm_dcm_x_neural(M.pE,model);
M.x   = x;
M.f   = f;
M.m   = size(M.pE.C,2);
M.g   = {};
M.ons = 32*ones(M.m,1);

% Volterra Kernels
%==========================================================================

% augment and bi-linearise
%--------------------------------------------------------------------------
[M0,M1,L1,L2] = spm_bireduce(M,M.pE);

% compute kernels (over 64 ms)
%--------------------------------------------------------------------------
M.ns    = 128;
U.dt    = pst/M.ns/1000;
t       = [1:M.ns]*U.dt*1000;
[K0,K1] = spm_kernels(M0,M1,L1,L2,M.ns,U.dt);

axes(handles.kernel)
plot(t,K1(:,:,1))
set(gca,'XLim',[min(t) max(t)])
ylabel('1st-order Volterra kernel')
xlabel('time (ms)')

% transfer function (source space)
%--------------------------------------------------------------------------
for i = 1:size(K1,2)
    g(:,i) = fft(K1(:,i));
    g(:,i) = abs(g(:,i).*conj(g(:,i)));
end
i   = [1:M.ns/2];
g   = g(i + 1,:);
w   = i/(M.ns*U.dt);

axes(handles.transfer)
plot(w,g)
set(gca,'XLim',[1 Hz])
xlabel('frequency time (Hz)')
ylabel('spectral density')

% prediction (source space)
%--------------------------------------------------------------------------
x     = spm_gen_erp(M.pE,M,U);
x     = x{1};

axes(handles.erp)
plot(t,x)
set(gca,'XLim',[min(t) max(t)])
xlabel('peristimulus time (ms)')
ylabel('voltage and conductance')

% fft (source space)
%--------------------------------------------------------------------------
clear g
for i = 1:size(x,2)
    g(:,i) = fft(x(:,i));
    g(:,i) = abs(g(:,i).*conj(g(:,i)));
end
i   = [1:M.ns/2];
g   = g(i + 1,:);
w   = i/(M.ns*U.dt);

axes(handles.fft)
plot(w,g)
set(gca,'XLim',[1 Hz])
xlabel('frequency time (Hz)')
ylabel('spectral density')

% unpack priors and display
%==========================================================================
function unpack_Callback(hObject, eventdata, handles);

% clear previous objects
%--------------------------------------------------------------------------
h = get(gcf,'Children');
for i = 1:length(h)
    if strcmp(get(h(i),'Tag'),'tmp')
        delete(h(i));
    end
end

% Type of model (neuronal)
%--------------------------------------------------------------------------
% 'ERP'    - (linear second order NMM slow)
% 'SEP'    - (linear second order NMM fast)
% 'CMC'    - (linear second order NMM Canonical microcircuit)
% 'LFP'    - (linear second order NMM self-inhibition)
% 'NMM'    - (nonlinear second order NMM first-order moments)
% 'MFM'    - (nonlinear second order NMM second-order moments)
try
    model = handles.DCM.options.model;
catch
    model = get(handles.model,'String');
    model = model{get(handles.model,'Value')};
    DCM.options.model = model;
end
switch model
    case{'ERP'}, set(handles.model,'Value',1);
    case{'SEP'}, set(handles.model,'Value',2);
    case{'CMC'}, set(handles.model,'Value',3);
    case{'LFP'}, set(handles.model,'Value',4);
    case{'NMM'}, set(handles.model,'Value',5);
    case{'MFM'}, set(handles.model,'Value',6);
    case{'NMDA'}, set(handles.model,'Value',8);
    otherwise
end

% get priors
%--------------------------------------------------------------------------
try
    handles.DCM.M.pE;
catch
    handles = reset_Callback(hObject, eventdata, handles);
end
pE = handles.DCM.M.pE;
pC = handles.DCM.M.pC;
try
    pC = spm_unvec(diag(pC),pE);
end

% display connection switches later
%--------------------------------------------------------------------------
try, pE = rmfield(pE,'B');,  end
try, pE = rmfield(pE,'SA');, end
try, pE = rmfield(pE,'GE');, end
try, pE = rmfield(pE,'GI');, end

% do not display spatial and spectral priors
%--------------------------------------------------------------------------
try, pE = rmfield(pE,{'Lpos','L','J'});,      end
try, pE = rmfield(pE,{'a','b','c','d'});, end

% display fields
%--------------------------------------------------------------------------
color = {[1 1 1],get(handles.name,'BackgroundColor')};

x0    = 1/8;
y0    = 1 - 1/8;
sx    = 1/24;
sy    = 1/48;
dx    = 1/256 + sx;
dy    = 1/256 + sy;
F     = fieldnames(pE);
for f = 1:length(F)

    P = getfield(pE,F{f});

    if iscell(P)

        % cell
        %------------------------------------------------------------------
        for k = 1:length(P)
            for i = 1:size(P{k},1);
                for j = 1:size(P{k},2)
                    x   = x0 + j*dx + (size(P{k},2) + 1)*dx*(k - 1);
                    y   = y0 - (i - 1)*dy;
                    str = sprintf('handles.DCM.M.pE.%s{%i}(%i,%i)',F{f},k,i,j);
                    str = ['handles=guidata(gcbo);' str '=str2num(get(gcbo,''String''));guidata(gcbo,handles);'];
                    h = uicontrol(gcf,...
                        'Units','Normalized',...
                        'Position',[x y sx sy],...
                        'Style','edit',...
                        'String',sprintf('%-4.2f',P{k}(i,j)),...
                        'enable','on',...
                        'Tag','tmp',...
                        'BackgroundColor',color{1 + ~eval(sprintf('pC.%s{%i}(%i,%i)',F{f},k,i,j))},...
                        'Callback',str);
                end
            end
        end

    elseif isvector(P)

        % vector
        %------------------------------------------------------------------
        for i   = 1:length(P);
            x   = x0 + i*dx;
            y   = y0;
            str = sprintf('handles.DCM.M.pE.%s(%i)',F{f},i);
            str = ['handles=guidata(gcbo);' str '=str2num(get(gcbo,''String''));guidata(gcbo,handles);'];
            uicontrol(gcf,...
                'Units','Normalized',...
                'Position',[x y sx sy],...
                'Style','edit',...
                'String',sprintf('%-4.2f',full(P(i))),...
                'enable','on',...
                'Tag','tmp',...
                'BackgroundColor',color{1 + ~eval(sprintf('pC.%s(%i)',F{f},i))},...
                'Callback',str);
        end
    else

        % matrix
        %------------------------------------------------------------------
        for i = 1:size(P,1);
            for j   = 1:size(P,2)
                x   = x0 + j*dx;
                y   = y0 - (i - 1)*dy;
                str = sprintf('handles.DCM.M.pE.%s(%i,%i)',F{f},i,j);
                str = ['handles=guidata(gcbo);' str '=str2num(get(gcbo,''String''));guidata(gcbo,handles);'];
                uicontrol(gcf,...
                    'Units','Normalized',...
                    'Position',[x y sx sy],...
                    'Style','edit',...
                    'String',sprintf('%-4.2f',full(P(i,j))),...
                    'enable','on',...
                    'Tag','tmp',...
                    'BackgroundColor',color{1 + ~eval(sprintf('pC.%s(%i,%i)',F{f},i,j))},...
                    'Callback',str);
            end
        end
    end

    % label
    %----------------------------------------------------------------------
    uicontrol(gcf,...
        'Units','Normalized',...
        'Position',[x0 - dx y0 2*sx sy],...
        'Style','text',...
        'String',sprintf('pE.%s',F{f}),...
        'ForegroundColor',[1 0 0],...
        'HorizontalAlignment','left',...
        'FontSize',12,...
        'Tag','tmp');
    y0 = y - dy - dy/2;

end

% connectivity matrices
%==========================================================================
pE      = handles.DCM.M.pE;
try, SA = pE.SA; catch, SA = []; end
try, GE = pE.GE; catch, GE = []; end
try, GI = pE.GI; catch, GI = []; end

% SA matrix
%--------------------------------------------------------------------------
x0    = 0.7;
y0    = 0.3;
for i = 1:size(SA,1);
    for j   = 1:size(SA,2)
        x   = x0 + j*dx;
        y   = y0 - (i - 1)*dy;
        str = sprintf('handles.DCM.M.pE.SA(%i,%i)',i,j);
        str = ['handles=guidata(gcbo);' str '=str2num(get(gcbo,''String''));guidata(gcbo,handles);'];
        uicontrol(gcf,...
            'Units','Normalized',...
            'Position',[x y sx sy],...
            'Style','edit',...
            'String',sprintf('%-4.2f',SA(i,j)),...
            'enable','on',...
            'Tag','tmp',...
            'Callback',str);
    end
end

% GE matrix
%--------------------------------------------------------------------------
x0    = 0.7;
y0    = 0.2;
for i = 1:size(GE,1);
    for j   = 1:size(GE,2)
        x   = x0 + j*dx;
        y   = y0 - (i - 1)*dy;
        str = sprintf('handles.DCM.M.pE.GE(%i,%i)',i,j);
        str = ['handles=guidata(gcbo);' str '=str2num(get(gcbo,''String''));guidata(gcbo,handles);'];
        uicontrol(gcf,...
            'Units','Normalized',...
            'Position',[x y sx sy],...
            'Style','edit',...
            'String',sprintf('%-4.2f',GE(i,j)),...
            'enable','on',...
            'Tag','tmp',...
            'Callback',str);
    end
end

% GI matrix
%--------------------------------------------------------------------------
x0    = 0.7;
y0    = 0.1;
for i = 1:size(GI,1);
    for j   = 1:size(GI,2)
        x   = x0 + j*dx;
        y   = y0 - (i - 1)*dy;
        str = sprintf('handles.DCM.M.pE.GI(%i,%i)',i,j);
        str = ['handles=guidata(gcbo);' str '=str2num(get(gcbo,''String''));guidata(gcbo,handles);'];
        uicontrol(gcf,...
            'Units','Normalized',...
            'Position',[x y sx sy],...
            'Style','edit',...
            'String',sprintf('%-4.2f',GI(i,j)),...
            'enable','on',...
            'Tag','tmp',...
            'Callback',str);
    end
end

% onset time
%==========================================================================
try
    ons  = handles.DCM.M.ons;
catch
    ons  = 32;
    handles.DCM.M.ons = ons;
end
str = 'handles.DCM.M.ons';
str = ['handles=guidata(gcbo);' str '=str2num(get(gcbo,''String''));guidata(gcbo,handles);'];
uicontrol(gcf,...
    'Units','Normalized',...
    'Position',[0.848 0.341 0.071 0.022],...
    'Style','edit',...
    'String',sprintf('%-4.0f',ons),...
    'enable','on',...
    'Tag','tmp',...
    'Callback',str);

guidata(hObject,handles);

% --- Executes on button press in reset.
%--------------------------------------------------------------------------
function handles = reset_Callback(hObject, eventdata, handles)

try
    DCM          = handles.DCM;
catch
    load_Callback(hObject, eventdata, handles)
    return
end
model            = DCM.options.model;
[pE,pC]          = spm_dcm_neural_priors(DCM.A,DCM.B,DCM.C,model);
handles.DCM.M.pE = pE;
handles.DCM.M.pC = pC;
guidata(hObject,handles);
unpack_Callback(hObject, eventdata, handles);

% --- Executes on button press in conditional.
%--------------------------------------------------------------------------
function conditional_Callback(hObject, eventdata, handles)

% conditional moments on parameters
%--------------------------------------------------------------------------
try
    handles.DCM.M.pE = handles.DCM.Ep;
    unpack_Callback(hObject, eventdata, handles);
catch
    warndlg('please invert and load a DCM to obtain conditional estimates')
end

% --- Executes on selection change in model.
function model_Callback(hObject, eventdata, handles)

% model type
%--------------------------------------------------------------------------
model = get(handles.model,       'String');
model = model{get(handles.model, 'Value')};
handles.DCM.options.model = model;
try
    reset_Callback(hObject, eventdata, handles);
catch
    load_Callback(hObject, eventdata, handles)
    return
end

author: spm


text:
function [D,P] = dijk(A,s,t)
%DIJK Shortest paths from nodes 's' to nodes 't' using Dijkstra algorithm.
% [D,p] = dijk(A,s,t)
%     A = n x n node-node weighted adjacency matrix of arc lengths
%         (Note: A(i,j) = 0   => Arc (i,j) does not exist;
%                A(i,j) = NaN => Arc (i,j) exists with 0 weight)
%     s = FROM node indices
%       = [] (default), paths from all nodes
%     t = TO node indices
%       = [] (default), paths to all nodes
%     D = |s| x |t| matrix of shortest path distances from 's' to 't'
%       = [D(i,j)], where D(i,j) = distance from node 'i' to node 'j'
%     P = |s| x n matrix of predecessor indices, where P(i,j) is the
%         index of the predecessor to node 'j' on the path from 's(i)' to 'j'
%         (use PRED2PATH to convert P to paths)
%       = path from 's' to 't', if |s| = |t| = 1
%
%  (If A is a triangular matrix, then computationally intensive node
%   selection step not needed since graph is acyclic (triangularity is a
%   sufficient, but not a necessary, condition for a graph to be acyclic)
%   and A can have non-negative elements)
%
%  (If |s| >> |t|, then DIJK is faster if DIJK(A',t,s) used, where D is now
%   transposed and P now represents successor indices)
%
%  (Based on Fig. 4.6 in Ahuja, Magnanti, and Orlin, Network Flows,
%   Prentice-Hall, 1993, p. 109.)

% Copyright (c) 1998-2001 by Michael G. Kay
% Matlog Version 5 22-Aug-2001

% Input Error Checking ******************************************************
error(nargchk(1,3,nargin));

[n,cA] = size(A);

if nargin < 2 | isempty(s), s = (1:n)'; else s = s(:); end
if nargin < 3 | isempty(t), t = (1:n)'; else t = t(:); end

if ~any(any(tril(A) ~= 0))       % A is upper triangular
   isAcyclic = 1;
elseif ~any(any(triu(A) ~= 0))   % A is lower triangular
   isAcyclic = 2;
else                             % Graph may not be acyclic
   isAcyclic = 0;
end

if n ~= cA
   error('A must be a square matrix');
elseif ~isAcyclic & any(any(A < 0))
   error('A must be non-negative');
elseif any(s < 1 | s > n)
   error(['''s'' must be an integer between 1 and ',num2str(n)]);
elseif any(t < 1 | t > n)
   error(['''t'' must be an integer between 1 and ',num2str(n)]);
end
% End (Input Error Checking) ************************************************

A = A';    % Use transpose to speed-up FIND for sparse A

D = zeros(length(s),length(t));
if nargout > 1, P = zeros(length(s),n); end

for i = 1:length(s)
   j = s(i);

   Di = Inf*ones(n,1); Di(j) = 0;

   isLab = logical(zeros(length(t),1));
   if isAcyclic ==  1
      nLab = j - 1;
   elseif isAcyclic == 2
      nLab = n - j;
   else
      nLab = 0;
      UnLab = 1:n;
      isUnLab = logical(ones(n,1));
   end

   while nLab < n & ~all(isLab)
      if isAcyclic
         Dj = Di(j);
      else	% Node selection
         [Dj,jj] = min(Di(isUnLab));
         j = UnLab(jj);
         UnLab(jj) = [];
         isUnLab(j) = 0;
      end

      nLab = nLab + 1;
      if length(t) < n, isLab = isLab | (j == t); end

      [jA,kA,Aj] = find(A(:,j));
      Aj(isnan(Aj)) = 0;

      if isempty(Aj), Dk = Inf; else Dk = Dj + Aj; end

      if nargout > 1, P(i,jA(Dk < Di(jA))) = j; end
      Di(jA) = min(Di(jA),Dk);

      if isAcyclic == 1       % Increment node index for upper triangular A
         j = j + 1;
      elseif isAcyclic == 2   % Decrement node index for lower triangular A
         j = j - 1;
      end
   end
   D(i,:) = Di(t)';
end

if nargout > 1 & length(s) == 1 & length(t) == 1
   P = pred2path(P,s,t);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% This file is part of the Matlab Toolbox for Dimensionality Reduction v0.7.2b.
% The toolbox can be obtained from http://homepage.tudelft.nl/19j49
% You are free to use, change, or redistribute this code in any way you
% want for non-commercial purposes. However, it is appreciated if you
% maintain the name of the original author.
%
% (C) Laurens van der Maaten, 2010
% University California, San Diego / Delft University of Technology

function scatter12n(name,data,labels)

ld=length(data(1,:));
    if ld==2
        hf=figure;
        set(hf,'name',name,'NumberTitle','off');
        %if handles.islb
            scatter(data(:,1),data(:,2),5,labels);
%         else
%             scatter(data(:,1),data(:,2),5);
%         end
        title(name);
    else
        if ld==1
            hf=figure;
            set(hf,'name',name,'NumberTitle','off');
            %if handles.islb
                scatter(data(:,1),zeros(length(data(:,1)),1),5,labels);
%             else
%                 scatter(data(:,1),zeros(length(data(:,1)),1),5);
%             end
            title(name);
        else
            %if handles.islb
                scattern(name,data,labels);
%             else
%                 scattern('Result of dimensionality reduction',data);
%             end
        end

    end

repo: cluster


text:
% DEMO_TEST.m

addpath(genpath('/home/vision/wren/caffe-dilate'));
caffe.reset_all();
clc; clear;

hazy_path = ('./inputs/');
type = '*.png';
hazy_data = dir(fullfile(hazy_path,type));

model_path = './models/';

solver_file = fullfile(model_path, 'dehaze_solver_test.prototxt');
save_file = fullfile(model_path, 'dehaze_ps128_bs1.mat');

Solver = modelconfig_test( solver_file, save_file);
% Solver.iter
Solver.batchsize = 1;
for ii = 45:length(hazy_data)
    disp(ii)
    clear batch;
    hazyimg = im2double(imread(fullfile(hazy_path,hazy_data(ii).name)));

    [row, col, cha] = size(hazyimg);

% resize the input to a multiple of 8
    height = ceil(row/8)*8;
    width = ceil(col/8)*8;
    hazyimg = imresize(hazyimg, [height, width]);

    hazy_wb = RealGWbal(uint8(255*hazyimg));
    hazy_wb = hazy_wb/255;
    hazy_cont = (2*(0.5+mean(hazyimg(:)))).*(hazyimg-mean(hazyimg(:)));
    hazy_gamma = hazyimg.^2.5;

    batch(:,:,1:3) = hazyimg;
    batch(:,:,4:6) = hazy_wb;
    batch(:,:,7:9) = hazy_cont;
    batch(:,:,10:12) = hazy_gamma;

%    figure(1);subplot(221), imshow(hazyimg);subplot(222),imshow(hazy_wb);
%              subplot(223), imshow(hazy_cont); subplot(224),imshow(hazy_gamma);

    batchc = {single(batch)};

    Solver.Solver_.net.blobs('data').reshape([height, width, 12, 1]);

    tic
    activec = Solver.Solver_.net.forward(batchc);
    toc
    dehazed= activec{3};

% 	dehazed = imresize(dehazed, [size(TEST_IMAGE, 1), size(TEST_IMAGE, 2)]);

	imwrite(dehazed, strcat('./results/',hazy_data(ii).name(1:end-4),'_dehazed.png'));
end
%imshow(out);title('Original/Proposed/GT');
%weights=Solver.Solver_.net.get_weights();
%data=Solver.Solver_.net.get_data();
%diff=Solver.Solver_.net.get_diff();

author: rwenqi

repo: GFN-dehazing

save_path: github-repos/MATLAB/rwenqi-GFN-dehazing


text:
function ensureCompatibleSymmetries(SO3VF,varargin)
% For calculating with @SO3VectorFields (+, -, .*, ./, ...) we have to verify
% that the symmetries are suitable.
%
% Syntax
%   ensureCompatibleSymmetries(SO3F1,SO3F2)
%   ensureCompatibleSymmetries(SO3F1,sF)
%   ensureCompatibleSymmetries(SO3F1,ori)
%   ensureCompatibleSymmetries(SO3F1,SO3F2,'conv')
%   ensureCompatibleSymmetries(SO3F1,'antipodal')
%
% Input
%  SO3F1, SO3F2 - @SO3VectorField, @SO3Fun
%  sF - @S2FunHarmonicSym
%  ori - @orientation
%
% Output
%  msg - yields a error message if the symmetry do not match
%
% Options
%  conv - be shure switched symmetries match
%

if isa(SO3VF,'SO3VectorField')
  SO3VF = SO3FunHarmonic(0,SO3VF.SRight,SO3VF.SLeft);
end
if nargin>1 && isa(varargin{1},'SO3VectorField')
  varargin{1} = SO3FunHarmonic(0,varargin{1}.SRight,varargin{1}.SLeft);
end

ensureCompatibleSymmetries(SO3VF,varargin{:},'SO3VectorField')

end

author: mtex-toolbox

repo: mtex

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


text:
function pen = sedumi2pen(F_struc,K,c,x0)
%SEDUMI2MPEN Internal function to convert SeDuMi structure to format needed in PENNON

% General data
pen.vars = length(c);
pen.constr = K.l;
if length(K.s)>1
    pen.mconstr = length(K.s);
else
    if K.s==0
        pen.mconstr=0;
    else
        pen.mconstr=1;
    end
end
pen.msizes = K.s;
pen.fobj = c(:)';
pen.x0 = zeros(1,length(c));

% Linear constraints
if K.l>0
    pen.ci = full(F_struc(K.f+1:K.f+K.l,1))';
    pen.bi_dim = zeros(1,K.f+K.l);
    [variables,constr,val] = find(-F_struc(1:K.l,2:end)');
    pen.bi_idx = variables(:)'-1;
    pen.bi_val = val(:)';
    j = 1;

    if ~isempty(constr)
        rr = histc(constr,[1:max(constr)]);
        uu = find(rr);
        pen.bi_dim(uu) = rr(uu);
        %uns = uniquestripped(constr);
        %for j = 1:length(uns)
        %    pen.bi_dim(uns(j)) = sum(constr==uns(j));
        %end
    end

else
    pen.ci = 0;
    pen.bi_dim = 0;
    pen.bi_idx = 0;
    pen.bi_val = 0;
end

% Semidefnite constraints
if K.s>0
    top = K.l+K.f+1;
    constraints = 1;
    pen.ai_dim = zeros(1,length(K.s));

    % First, optimized code for scalar SDPs (happens for nonlinear scalar constraints)
    if all(K.s==1)

        Avec = -F_struc(top:top+length(K.s)-1,:);
        [variables,constr,val] = find(Avec');
        j = 1;
        n_nz = length(val);
        pen.ai_nzs = repmat(1,1,n_nz);
        pen.ai_row = repmat(0,1,n_nz);
        pen.ai_col = pen.ai_row;%repmat(0,1,n_nz);
        pen.ai_val = val(:)';
        pen.ai_idx = variables(:)'-1;
        while j<=length(val)
            first_constraint = constr(j);
            n_nz = 1;
            while j+n_nz<=length(val) & constr(j+n_nz)==first_constraint
                n_nz = n_nz + 1;
            end
            pen.ai_dim(first_constraint) = n_nz;
            j = j + n_nz;
        end
    else
        pen.ai_idx = [];
        pen.ai_val = [];
        pen.ai_row = [];
        pen.ai_col = [];
        pen.ai_nzs = [];
        while (constraints<=length(K.s))
            n = K.s(constraints);
            Avec = -F_struc(top:top+n^2-1,:);
            picker = reshape(1:n^2,n,n);picker = tril(picker-diag(diag(picker)));picker = find(picker(:));
            [rowcols,varindicies,vals]=find(Avec);
            %the_col = 1+floor((rowcols-1)/n);
            %the_row = rowcols-(the_col-1)*n;
            %removethese = the_row > the_col;
            removethese = find(ismembc(rowcols,picker));
            rowcols(removethese) =[];
            varindicies(removethese)=[];
            vals(removethese)=[];
            if ~isempty(rowcols)
                cols = ceil(rowcols/n);
                rows = rowcols - n*(cols-1);
                uns = uniquestripped(varindicies);

                difference = diff([varindicies(:) ; max(varindicies)+1]);
                count = diff(find([1 ; difference]))';
                pen.ai_nzs = [pen.ai_nzs count];

                pen.ai_dim(constraints) =  length(uniquestripped(varindicies));
                pen.ai_idx = [pen.ai_idx uns(:)'-1];
                pen.ai_row = [pen.ai_row rows(:)'-1];
                pen.ai_col = [pen.ai_col cols(:)'-1];
                pen.ai_val = [pen.ai_val vals(:)'];
            end
            constraints = constraints+1;
            top = top+n*n;
        end
    end
else
    pen.ai_dim = 0;
    pen.ai_idx = 0;
    pen.ai_val = 0;
    pen.ai_row = 0;
    pen.ai_col = 0;
    pen.ai_nzs = 0;
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function euler_number2_test ( )

%*****************************************************************************80
%
%% POLPAK_TEST0265 tests EULER_NUMBER2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'EULER_NUMBER2_TEST\n' );
  fprintf ( 1, '  EULER_NUMBER2 computes Euler numbers.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  N  exact   EULER_NUMBER2\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, c1 ] = euler_number_values ( n_data );

    if ( n_data == 0 )
      break
    end

    c2 = euler_number2 ( n );

    fprintf ( 1, '  %2d  %10d  %12f\n', n, c1, c2 );

  end

  return
end


text:
function x = meyer(p, data1, data2)
  n=16;

% data1, data2 are actually unused

  for i=1:n
    ui=0.45+0.05*i;
    x(i)=p(1)*exp(10.0*p(2)/(ui+p(3)) - 13.0);
  end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
RadarName = 'staddon_detection';
RadarCoords.lat = 50.346069;
RadarCoords.lon = -4.113670;

RadarCoords_lat = 50.339026;
RadarCoords_lon = -4.126307;

% Load dataset
load(strcat(RadarName,'.mat'));
% load('jpda_mon_no_detection.mat');

% lon_lim = [-4.2637   -4.1024];
% lat_lim = [50.2699   50.3853];

lon_lim_obs = [-4.16737200000000,-4.14897200000000];
lat_lim_obs = [50.3537860000000,50.3617750000000];
% radar_meas_convert;

% Map plot
figure('units','normalized','outerposition',[0 0 .5 0.8])
ax(1) = gca;
plot_google_map('Axis',ax(1),'APIKey','AIzaSyBXKujdtXRZiqya1soVS9pxBzYR4g7aGvM','Resize',3,'Scale',2,'MapType','satellite');
axis(ax(1),[lon_lim(1) lon_lim(2) lat_lim(1) lat_lim(2)])
plots = [];
hold on;

N = size(MeasurementScans,2); % Simulation length

Tracks = [TrackList, DeletedTracks];
t = [368, 912, 640, 673, 146, 594, 342, 182, 912, 505, 672, 195, 594, 342, ...
    798, 267, 54, 794, 826, 125, 767, 524, 63, 200, 754, 244, 74, 777, 969, ...
    77, 51, 74, 828, 307, 861, 97, 807, 931, 665, 404, 177, 938, 225, 392, ...
    985, 585, 608, 319, 641, 915, 661, 538, 517, 586, 27, 164];
for k=2:N

    MeasurementList = MeasurementScans(k);
    if(MeasurementList.NumMeasurements>0)
        % Convert measurements to LLA and plot them
    %     [lat,lon,~] = ned2geodetic(MeasurementList.Vectors(2,:),...
    %                                MeasurementList.Vectors(1,:),...
    %                                0,...
    %                                RadarCoords.lat,...
    %                                RadarCoords.lon,...
    %                                0,...
    %                                referenceEllipsoid('wgs84'));
        latlon = [MeasurementList.Metadata.LatLon];
        lat = latlon(1,:);
        lon = latlon(2,:);
        plots(end+1) = plot(ax(1), lon,lat,'y.','MarkerSize', 4);
    end
end

% Plot all existing tracks
for j=1:numel(Tracks)
    track = Tracks{j};
%     if find(t==track.Tag.ID,1)
%         continue
%     end
    % Convert track trajectory to LLA and plot it
    means = [track.Trajectory.Mean];
    [lat,lon,~] = ned2geodetic(means(3,:),...
                               means(1,:),...
                               0,...
                               RadarCoords.lat,...
                               RadarCoords.lon,...
                               0,...
                               referenceEllipsoid('wgs84'));
    traj_length = size(lon,2);
    if(traj_length>NumPersistFrames)
        start = traj_length-NumPersistFrames;
    else
        start = 1;
    end
    plots(end+1) = plot(ax(1), lon(:,start),lat(:,start),'go','MarkerSize',15,'LineWidth',2);
    plots(end+1) = plot(ax(1), lon(:,start:end),lat(:,start:end),'-.w','LineWidth',2);
    plots(end+1) = plot(ax(1), lon(:,end),lat(:,end),'rx','MarkerSize',15,'LineWidth',2);

    % Convert track velocity to LLA and plot it
%     [lat_vel,lon_vel,~] = ned2geodetic(track.State.Mean(4,end),...
%                                        track.State.Mean(2,end),...
%                                        0,...
%                                        lat(:,end),...
%                                        lon(:,end),...
%                                        0,...
%                                        referenceEllipsoid('wgs84'));
%     lat_vel = lat_vel-lat(:,end);
%     lon_vel = lon_vel-lon(:,end);
%     plots(end+1) = quiver(ax(1), lon(:,end),lat(:,end),20*lon_vel,20*lat_vel,'r','LineWidth',1.5);
%     if ShowTrackInfo
%     text(ax(1), lon(:,end), lat(:,end), num2str(track.Tag.ID),'FontSize',18,'Color','g')
%     end

end

% Plot region of low PD
x = [-3422.85261310741,-3191.58441229017,-2993.55608122595,-3228.08558459284, -3422.85261310741];
y = [1.472966334316646e+03,1.123217741935625e+03,1.243952901078573e+03,1.579667558371468e+03, 1.472966334316646e+03];
[y,x,~] = ned2geodetic(y,...
                       x,...
                       0,...
                       RadarCoords.lat,...
                       RadarCoords.lon,...
                       0,...
                       referenceEllipsoid('wgs84'));
plots(end+1) = plot(ax(1), x, y, 'm-', 'LineWidth', 2);

% Plot region of high clutter
x = [-469.710602736631,-2572.92295458513,-2516.31234127506,-1320.74163107899,-835.307363807233, -469.710602736631];
y = [-4915.41679013513,-4065.50122858976,-2858.60918927224,-1380.19032924249,-1261.08679763585, -4915.41679013513];
%         x = [-2858.62556879160,-2856.97056906902,-108.925126527319,-108.988225201227,-2858.62556879160];
%         y = [-4044.74840882411,-974.655039469575,-975.424226884065,-4045.51804181679,-4044.74840882411];
[y,x,~] = ned2geodetic(y,...
                       x,...
                       0,...
                       RadarCoords.lat,...
                       RadarCoords.lon,...
                       0,...
                       referenceEllipsoid('wgs84'));
plots(end+1) = plot(ax(1), x, y, 'r-', 'LineWidth', 2);

plot(ax(1), RadarCoords_lon,RadarCoords_lat,...
                               '-s','MarkerSize',10,...
                               'MarkerEdgeColor','red',...
                               'MarkerFaceColor',[1 .6 .6]);

h1 = plot(NaN, NaN, 'o','MarkerSize', 5, 'MarkerFaceColor', 'y', 'MarkerEdgeColor','k');
h2 = plot(NaN, NaN, 's','MarkerSize',10,...
                               'MarkerEdgeColor','red',...
                               'MarkerFaceColor',[1 .6 .6]);
h3 = plot(NaN, NaN, 'm-', 'LineWidth', 2);
h4 = plot(NaN, NaN, 'r-', 'LineWidth', 2);
h5 = plot(NaN, NaN, 'go','MarkerSize',15,'LineWidth',2);
h6 = plot(NaN, NaN, 'rx','MarkerSize',15,'LineWidth',2);
legend([h1, h2, h5, h6], 'Detections', 'Obscured Region', 'Track Birth', 'Track Death');
% legend([h1, h2, h3, h4], 'Detections', 'Radar Position', 'Obscured Region', 'High-clutter region');
% xlabel('Longitude');
% ylabel('Latitude');
set(ax(1),'YTickLabel',[]);
set(ax(1),'XTickLabel',[]);

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function cS = antigorite
% serpentine

cs_Atg = crystalSymmetry('m',[4.7013 1 0.7844], [90, 91.633, 90]*degree);
N = Miller({0,0,1},{0,0,-1},{1,0,0},{-1,0,0},{2,1,0},{-2,1,0},cs_Atg);
dist = [0.15, 0.15, 0.4, 0.4, 2, 2];
cS = crystalShape(N./dist);


text:
function oneDfeature = trans2Dto1DFunc(twoDfeature)
rowSize = size(twoDfeature,1);
oneDfeature = [];
for i=1:rowSize
    oneDfeature = [oneDfeature,twoDfeature(i,:)];
end

author: DrGabor

repo: LiDAR

save_path: github-repos/MATLAB/DrGabor-LiDAR


text:
function varargout = delLinks(varargin)
% DELLINKS M-file for delLinks.fig
%      DELLINKS, by itself, creates a new DELLINKS or raises the existing
%      singleton*.
%
%      H = DELLINKS returns the handle to a new DELLINKS or the handle to
%      the existing singleton*.
%
%      DELLINKS('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in DELLINKS.M with the given input arguments.
%
%      DELLINKS('Property','Value',...) creates a new DELLINKS or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before delLinks_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to delLinks_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help delLinks

% Last Modified by GUIDE v2.5 26-Oct-2006 13:44:29

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @delLinks_OpeningFcn, ...
                   'gui_OutputFcn',  @delLinks_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin & isstr(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before delLinks is made visible.
function delLinks_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to delLinks (see VARARGIN)

% Choose default command line output for delLinks
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes delLinks wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = delLinks_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes during object creation, after setting all properties.
function EdSelectedSubsystem_CreateFcn(hObject, eventdata, handles)
% hObject    handle to EdSelectedSubsystem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function EdSelectedSubsystem_Callback(hObject, eventdata, handles)
% hObject    handle to EdSelectedSubsystem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of EdSelectedSubsystem as text
%        str2double(get(hObject,'String')) returns contents of EdSelectedSubsystem as a double

% --- Executes on button press in PbSelectSubsystem.
function PbSelectSubsystem_Callback(hObject, eventdata, handles)
% hObject    handle to PbSelectSubsystem (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
set(handles.EdSelectedSubsystem,'String', gcb);

% --- Executes on button press in PbDeleteLinks.
function PbDeleteLinks_Callback(hObject, eventdata, handles)
% hObject    handle to PbDeleteLinks (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

NumDelLinks=0;

if (strcmp(get(handles.EdSelectedSubsystem,'String'),'No Selection')==1)
    errordlg('Select a Subsystem first!','No Subsystem Selected');
    return
else
    answer = questdlg('All links in the selected subsystem will be broken! Continue?',...
        get(handles.EdSelectedSubsystem,'String'),'Yes','No','No');

    if strcmp (answer,'Yes')
        libs={'dummy'};

        while isempty(libs)==0
            libs=find_system(gcb,'LookUnderMasks','all','LinkStatus','resolved');
            libStruct=cell2struct(libs,'libPath',2);
            cntLibPath=size(libStruct);

            for i=1:cntLibPath(1)
                set_param(libStruct(i).libPath,'LinkStatus','none');
                NumDelLinks=NumDelLinks+1;
            end
        end
        warndlg([num2str(NumDelLinks) ,' link(s) broken. Changes take effect after you have saved the model!'],'Links broken');
        return;
    else
        return ;
    end;
end;
%if needed:
%LinkStatus 'resolved' = Library Link active
%LinkStatus 'inactive' = Library Link disabled
%LinkStatus 'none' = Library Link broken

author: Sable

repo: mcbench-benchmarks


text:
function [A,b,x] = shaw(n)
%SHAW Test problem: one-dimensional image restoration model.
%
% [A,b,x] = shaw(n)
%
% Discretization of a first kind Fredholm integral equation with
% [-pi/2,pi/2] as both integration intervals.  The kernel K and
% the solution f, which are given by
%    K(s,t) = (cos(s) + cos(t))*(sin(u)/u)^2
%    u = pi*(sin(s) + sin(t))
%    f(t) = a1*exp(-c1*(t - t1)^2) + a2*exp(-c2*(t - t2)^2) ,
% are discretized by simple quadrature to produce A and x.
% Then the discrete right-hand b side is produced as b = A*x.
%
% The order n must be even.

% Reference: C. B. Shaw, Jr., "Improvements of the resolution of
% an instrument by numerical solution of an integral equation",
% J. Math. Anal. Appl. 37 (1972), 83-112.

% Per Christian Hansen, IMM, 08/20/91.

% Check input.
if (rem(n,2)~=0), error('The order n must be even'), end

% Initialization.
h = pi/n; A = zeros(n,n);

% Compute the matrix A.
co = cos(-pi/2 + (.5:n-.5)*h);
psi = pi*sin(-pi/2 + (.5:n-.5)*h);
for i=1:n/2
  for j=i:n-i
    ss = psi(i) + psi(j);
    A(i,j) = ((co(i) + co(j))*sin(ss)/ss)^2;
    A(n-j+1,n-i+1) = A(i,j);
  end
  A(i,n-i+1) = (2*co(i))^2;
end
A = A + triu(A,1)'; A = A*h;

% Compute the vectors x and b.
a1 = 2; c1 = 6; t1 =  .8;
a2 = 1; c2 = 2; t2 = -.5;
if (nargout>1)
  x =   a1*exp(-c1*(-pi/2 + (.5:n-.5)'*h - t1).^2) ...
      + a2*exp(-c2*(-pi/2 + (.5:n-.5)'*h - t2).^2);
  b = A*x;
end

repo: REST


text:
function [mus,sigmas,R] = restore_from_projection(mus_p,sigmas_p,R_p,c,E)
% restore from projection.
% input:
%   mus_p - 1 by M cell array (each cell is a K_j by B matrix)
%   sigmas_p - 1 by M cell array (each cell is a B by B by K_j matrix)
%   R_p - M by d projected data or M by d by N projected data
%   c - 1 by B center of the PCA
%   E - B by d projection matrix
% output:
%   mus, sigmas, R - restored versions of mus_p, sigmas_p and R_p

if ~isempty(mus_p) && ~isempty(sigmas_p)
    M = length(mus_p);
elseif ~isempty(R_p)
    M = size(R_p,1);
end
mus = cell(1,M);
sigmas = cell(1,M);
B = size(E,1);
R = zeros(M,B);

if ~isempty(mus_p) && ~isempty(sigmas_p)
    for j = 1:M
        K_j = size(mus_p{j},1);
        mus{j} = zeros(K_j,B);
        sigmas{j} = zeros(B,B,K_j);

        for k = 1:size(mus_p{j},1)
    %         mus{j}(k,:) = solve_linsys_mu(E,mus_p{j}(k,:)')' + c;
    %         sigmas{j}(:,:,k) = solve_linsys_sigma(E, sigmas_p{j}(:,:,k));
            mus{j}(k,:) = recover_mu(E,c,mus_p{j}(k,:));
            sigmas{j}(:,:,k) = recover_sigma(E,sigmas_p{j}(:,:,k));
        end
    end
end

if ~isempty(R_p)
    if ndims(R_p) == 2
        R = recover_mu(E,c,R_p);
    elseif ndims(R_p) == 3
        R = zeros(M,B,size(R_p,3));
        for i = 1:size(R_p,3)
            R(:,:,i) = recover_mu(E,c,R_p(:,:,i));
        end
    end
end

function mu = recover_mu(E,c,mu0)
mu = E*mu0' + repmat(c',[1, size(mu0,1)]);
mu = mu';

function sigma = recover_sigma(E, sigma0)
sigma = E*sigma0*E' + 1e-9*eye(size(E,1));

function x = solve_linsys_mu(E,mu0)
EE1 = E*E';
B = size(EE1,1);
EE1 = EE1 + 1e-9*eye(B);
x = EE1 \ (E* mu0);

function sigma = solve_linsys_sigma(E,sigma0)
[B,d] = size(E);
w = 1e-9;
EE = kron(E,E);
tmp1 = EE * sigma0(:);
tmp2 = (eye(d^2) + (1/w)*EE'*EE) \ (EE' * tmp1);
sigma = (1/w) * tmp1 - (1/w^2) * EE * tmp2;
sigma = reshape(sigma,B,B);

author: zhouyuanzxcv

repo: Hyperspectral