text:
%% Copyright (C) 2014, 2016 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @defmethod  @@sym private_disp_name (@var{x}, @var{name})
%% A string appropriate for representing the name of this sym.
%%
%% Private method: this is not the method you are looking for.
%%
%% @end defmethod

function s = private_disp_name(x, input_name)

  s = input_name;
  % subclasses might do something more interesting, but they should
  % be careful to ensure empty input_name gives empty s.

end

%!test
%! syms x
%! s = private_disp_name(x, 'x');
%! assert (strcmp (s, 'x'))

repo: octsympy


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% thermalLB.m: Rayleigh Benard Convection, using a LB method,
%   based on [Z.Guo, e.a., http://dx.doi.org/10.1002/fld.337].
%   Boussinesq approximation is used for the buoyancy term:
%     - Fluid is approximated with incompressible Navier-Stokes
%       equations including a body force term, and simulated
%       with a BGK model
%     - Temperature is approximated with advection-diffusion
%       equation and simulated with a BGK model
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lattice Boltzmann sample, written in matlab
% Copyright (C) 2008 Andrea Parmigiani, Orestis Malaspinas, Jonas Latt
% Address: Rue General Dufour 24,  1211 Geneva 4, Switzerland
% E-mail: andrea.parmigiani@terre.unige.ch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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., 51 Franklin Street, Fifth Floor,
% Boston, MA  02110-1301, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear

% GENERAL FLOW CONSTANTS

ly           = 51;
aspect_ratio = 2;
lx           = aspect_ratio*ly;
delta_x      = 1./(ly-2);
Pr           = 1.;
Ra           = 20000.; % Rayleigh number
gr           = 0.001;  % Gravity
buoyancy     = [0,gr];

Thot  = 1; % Heating on bottom wall
Tcold = 0; % Cooling on top wall
T0 = (Thot+Tcold)/2;

delta_t = sqrt(gr*delta_x);
% nu: kinematic viscosity in lattice units
nu      = sqrt(Pr/Ra)*delta_t/(delta_x*delta_x);
% k: thermal diffusivity
k       = sqrt(1./(Pr*Ra))*delta_t/(delta_x*delta_x);
omegaNS = 1./(3*nu+0.5); % Relaxation parameter for fluid
omegaT  = 1./(3.*k+0.5); % Relaxation parameter for temperature

maxT   = 80000;    % total number of iterations
tPlot  = 100;      % iterations between successive graphical outputs
tStatistics = 10;  % iterations between successive file accesses

% D2Q9 LATTICE CONSTANTS
tNS   = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36];
cxNS  = [  0,   1,  0, -1,  0,    1,  -1,  -1,   1];
cyNS  = [  0,   0,  1,  0, -1,    1,   1,  -1,  -1];
oppNS = [  1,   4,  5,  2,  3,    8,   9,   6,   7];

% D2Q5 LATTICE CONSTANTS
tT   = [1/3, 1/6, 1/6, 1/6, 1/6];
cxT  = [  0,   1,   0,  -1,   0];
cyT  = [  0,   0,   1,   0,  -1];
oppT = [  1,   4,   5,   2,   3];

[y,x] = meshgrid(1:ly,1:lx);

% INITIAL CONDITION FOR FLUID: (rho=1, u=0) ==> fIn(i) = t(i)
fIn = reshape( tNS' * ones(1,lx*ly), 9, lx, ly);

% INITIAL CONDITION FOR TEMPERATURE: (T=0) ==> TIn(i) = t(i)
tIn = reshape( tT' *Tcold *ones(1,lx*ly), 5, lx, ly);
% Except for bottom wall, where T=1
tIn(:,:,ly)=Thot*tT'*ones(1,lx);
% We need a small trigger, to break symmetry
tIn(:,lx/2,ly-1)= tT*(Thot + (Thot/10.));

% Open file for statistics
fid = fopen('thermal_statistics.dat','w');
fprintf(fid,'Thermal Statistics: time-step --- uy[nx/2,ny/2] --- Nu\n\n\n');

% MAIN LOOP (TIME CYCLES)
for cycle = 1:maxT
  % MACROSCOPIC VARIABLES
   rho = sum(fIn);
   T = sum(tIn); %temperature
   ux  = reshape ( (cxNS * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho;
   uy  = reshape ( (cyNS * reshape(fIn,9,lx*ly)), 1,lx,ly) ./rho;

   % MACROSCOPIC BOUNDARY CONDITIONS
   % NO-SLIP for fluid and CONSTANT at lower and upper
   % boundary...  periodicity wrt. left-right
   % COLLISION STEP FLUID
   for i=1:9
      cuNS         = 3*(cxNS(i)*ux+cyNS(i)*uy);
      fEq(i,:,:)   = rho .* tNS(i) .* ...
                       ( 1 + cuNS + 1/2*(cuNS.*cuNS) - 3/2*(ux.^2+uy.^2) );
      force(i,:,:) = 3.*tNS(i) .*rho .* (T-T0) .* ...
                       (cxNS(i)*buoyancy(1)+cyNS(i)*buoyancy(2))/(Thot-Tcold);
      fOut(i,:,:)  = fIn(i,:,:) - omegaNS .* (fIn(i,:,:)-fEq(i,:,:)) + force(i,:,:);
   end

    % COLLISION STEP TEMPERATURE
   for i=1:5
      cu          = 3*(cxT(i)*ux+cyT(i)*uy);
      tEq(i,:,:)  = T .* tT(i) .* ( 1 + cu );
      tOut(i,:,:) = tIn(i,:,:) - omegaT .* (tIn(i,:,:)-tEq(i,:,:));
   end

   % MICROSCOPIC BOUNDARY CONDITIONS FOR FLUID
   for i=1:9
        fOut(i,:,1)  = fIn(oppNS(i),:,1);
        fOut(i,:,ly) = fIn(oppNS(i),:,ly);
   end

   % STREAMING STEP FLUID
   for i=1:9
      fIn(i,:,:) = circshift(fOut(i,:,:), [0,cxNS(i),cyNS(i)]);
   end

   % STREAMING STEP FLUID
   for i=1:5
      tIn(i,:,:) = circshift(tOut(i,:,:), [0,cxT(i),cyT(i)]);
   end

   % MICROSCOPIC BOUNDARY CONDITIONS FOR TEMEPERATURE
   %
   tIn(5,:,ly) = Tcold-tIn(1,:,ly)-tIn(2,:,ly)-tIn(3,:,ly)-tIn(4,:,ly);
   tIn(3,:,1)  = Thot-tIn(1,:,1)  -tIn(2,:,1) -tIn(4,:,1) -tIn(5,:,1);

   % VISUALIZATION
   if (mod(cycle,tStatistics)==0)
       u     = reshape(sqrt(ux.^2+uy.^2),lx,ly);
       uy_Nu = reshape(uy,lx,ly); % vertical velocity
       T     = reshape(T,lx,ly);
       Nu    = 1. + sum(sum(uy_Nu.*T))/(lx*k*(Thot-Tcold));
       fprintf(fid,'%8.0f  %12.8f  %12.8f\n',cycle,u(int8(lx/2),int8(ly/2))^2, Nu);
       if(mod(cycle,tPlot)==0)
           subplot(2,1,1);
           imagesc(u(:,ly:-1:1)');
           title('Fluid velocity');
           axis off; drawnow
           subplot(2,1,2);
           imagesc(T(:,ly:-1:1)')
           title(['Temperature (Nusselt number is ' num2str(Nu) ')']);
           axis off; drawnow
       end
   end
end

fclose(fid);

repo: Aero-matlab


text:
function [f,df] = sqplab_fun(x)

global SQPLABDATA

xevaled = zeros(1,length(SQPLABDATA.interfacedata.c));
xevaled(SQPLABDATA.linearindicies) = x;

% Experimental support for arbitrary functions
if SQPLABDATA.evalinobjective
    if SQPLABDATA.monominobjective
        if ~isempty(SQPLABDATA.bilinears)
            xevaled(SQPLABDATA.bilinears(:,1)) = xevaled(SQPLABDATA.bilinears(:,2)).*xevaled(SQPLABDATA.bilinears(:,3));
        else
            xevaled(SQPLABDATA.nonlinearindicies) = prod(repmat(xevaled,length(SQPLABDATA.nonlinearindicies),1).^SQPLABDATA.monomtable(SQPLABDATA.nonlinearindicies,:),2);
        end
    end
    for i = 1:length(SQPLABDATA.interfacedata.evalMap)
        arguments = SQPLABDATA.evalMap{i}.prearg;
        arguments{2} = xevaled(SQPLABDATA.interfacedata.evalMap{i}.variableIndex);
        xevaled(SQPLABDATA.interfacedata.evalVariables(i)) = feval(arguments{:});
    end
end

if SQPLABDATA.monominobjective
    if ~isempty(SQPLABDATA.bilinears)
        xevaled(SQPLABDATA.bilinears(:,1)) = xevaled(SQPLABDATA.bilinears(:,2)).*xevaled(SQPLABDATA.bilinears(:,3));
    else
        xevaled(SQPLABDATA.nonlinearindicies) = prod(repmat(xevaled,length(SQPLABDATA.nonlinearindicies),1).^SQPLABDATA.monomtable(SQPLABDATA.nonlinearindicies,:),2);
    end
end

xevaled = xevaled(:);
if SQPLABDATA.SimpleLinearObjective
    f = SQPLABDATA.interfacedata.c'*xevaled;
else
    f = (SQPLABDATA.interfacedata.c'+xevaled'*SQPLABDATA.interfacedata.Q)*xevaled;
end

if SQPLABDATA.SimpleLinearObjective
    df = SQPLABDATA.interfacedata.c(SQPLABDATA.linearindicies);
elseif SQPLABDATA.SimpleQuadraticObjective
    df = SQPLABDATA.interfacedata.c(SQPLABDATA.linearindicies) + 2*SQPLABDATA.interfacedata.Q(SQPLABDATA.linearindicies,SQPLABDATA.linearindicies)*x;
elseif SQPLABDATA.SimpleNonlinearObjective
    df = [];
    for i = 1:length(SQPLABDATA.linearindicies)
        xevaled = zeros(1,length(SQPLABDATA.interfacedata.c));
        xevaled(SQPLABDATA.linearindicies) = x;
        mt = SQPLABDATA.monomtable;
        oldpower = mt(:,SQPLABDATA.linearindicies(i));
        mt(:,SQPLABDATA.linearindicies(i)) = mt(:,SQPLABDATA.linearindicies(i))-1;
        xevaled = prod(repmat(xevaled,size(mt,1),1).^mt,2);
        xevaled = xevaled(:)'.*oldpower';xevaled(isnan(xevaled))=0;
        df = [df;SQPLABDATA.interfacedata.c'*xevaled'];
    end
    df = df + 2*SQPLABDATA.interfacedata.Q(SQPLABDATA.linearindicies,SQPLABDATA.linearindicies)*x;
else
    df = [];
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% SYNTAX:
% dod = hmrR_Intensity2OD( intensity )
%
% UI NAME:
% Intensity_to_Delta_OD
%
% DESCRIPTION:
% Converts intensity data to optical density
%
% INPUT:
% intensity - SNIRF data type where the d matrix is intensity
%
% OUTPUT:
% dod - SNIRF data type where the d matrix is change in optical density
%
% USAGE OPTIONS:
% Intensity_to_Delta_OD: dod = hmrR_Intensity2OD(data)
%
function dod = hmrR_Intensity2OD( intensity )

% convert to dod
dod = DataClass().empty();
for ii=1:length(intensity)
    dod(ii) = DataClass();
    d = intensity(ii).GetDataTimeSeries();
    dm = mean(abs(d),1);
    nTpts = size(d,1);
    dod(ii).SetTime(intensity(ii).GetTime());
    dod(ii).SetDataTimeSeries(-log(abs(d)./(ones(nTpts,1)*dm)));
    dod(ii).SetMl(intensity(ii).GetMl());
    dod(ii).SetDataTypeDod();
end

author: BUNPC


text:
function [y] = spm_gx_fmri_linear(x,u,P,M)
% Simulated BOLD response to input (linear version)
% FORMAT [y] = spm_gx_fmri_linear(x,u,P,M)
% y          - BOLD response (%)
% x          - state vector     (see spm_fx_fmri)
% P          - Parameter vector (see spm_fx_fmri)
% M          - model specification structure (see spm_nlsi)
%__________________________________________________________________________
%
% This function implements the BOLD signal model described in:
%
% Stephan KE, Weiskopf N, Drysdale PM, Robinson PA, Friston KJ (2007)
% Comparing hemodynamic models with DCM. NeuroImage 38: 387-401.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston & Klaas Enno Stephan
% $Id: spm_gx_fmri_linear.m 6262 2014-11-17 13:47:56Z karl $

% Biophysical constants for 1.5T
%==========================================================================

% time to echo (TE) (default 0.04 sec)
%--------------------------------------------------------------------------
try, TE = M.TE; catch, TE = 0.04; end

% resting venous volume (%)
%--------------------------------------------------------------------------
V0  = 4;

% estimated region-specific ratios of intra- to extra-vascular signal
%--------------------------------------------------------------------------
ep  = 1*exp(P.epsilon);

% slope r0 of intravascular relaxation rate R_iv as a function of oxygen
% saturation S:  R_iv = r0*[(1 - S)-(1 - S0)] (Hz)
%--------------------------------------------------------------------------
r0  = 25;

% frequency offset at the outer surface of magnetized vessels (Hz)
%--------------------------------------------------------------------------
nu0 = 40.3;

% resting oxygen extraction fraction
%--------------------------------------------------------------------------
E0  = 0.4;

%-Coefficients in BOLD signal model
%==========================================================================
k1  = 4.3*nu0*E0*TE;
k2  = ep*r0*E0*TE;
k3  = 1 - ep;

%-Output equation of BOLD signal model
%==========================================================================
v   = x(:,4) + 1;
q   = x(:,5) + 1;
y   = V0*((k1 + k2).*(1 - q) + (k3 - k2).*(1 - v));

author: spm


text:
function plotmimac(mi,inde)

    % TODO maybe move into domisfit

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun(mfilename('fullpath'));

    global mif2 mif1
    global oneOfHowManyPopupIdx

    %var1 = 4;
    sc = oneOfHowManyPopupIdx;
    figure(UNK) % FIXME: really? this figure? unsure
    delete(findobj(UNK,'Type','axes'));
    rect = [0.15,  0.20, 0.75, 0.65];
    axes('position',rect)
    watchon

    % check if cross-section exists
    figNumber=findobj('Type','Figure','-and','Name','Cross -Section');

    if isempty(figNumber)
        errordlg('Please create a cross-section first, then rerun the last selection');
        nlammap
        return
    end

    % check if cross-section is still current
    if max(mi(:,1)) > length(mi(:,1))
        errordlg('Please rerun the cross-section first, then rerun the last selection');
        nlammap
        return
    end

    mic = mi(inde,:);
    le = size(newa,2); %FIXME where does newa come from? ZG.newa? input parameter?  Needs to be treated like a ZmapCatalog

    if var1 == 1
        for i = 1:length(newa(:,6))
            pl =  plot(newa(i,le),-newa(i,7),'ro');
            hold on
            set(pl,'MarkerSize',mic(i,2)/sc)
        end

    elseif var1 == 2

        for i = 1:length(newa(:,6))
            pl =  plot(newa(i,le),-newa(i,7),'bx');
            hold on
            set(pl,'MarkerSize',mic(i,2)/sc,'LineWidth',mic(i,2)/sc)
        end

    elseif var1 == 3

        for i = 1:length(newa(:,6))
            pl =  plot(newa(i,le),-newa(i,7),'bx');
            hold on
            c = mic(i,2)/max(mic(:,2));
            %c = newa(i,15)*10;
            set(pl,'MarkerSize',mic(i,2)/sc+3,'LineWidth',mic(i,2)/sc+0.5,'Color',[ c c c ] )
        end

    elseif var1 == 4

        g = jet;
        for i = 1:length(newa(:,6))
            pl =  plot(newa(i,le),-newa(i,7),'bx');
            hold on
            c = floor(mic(i,2)/max(mic(:,2))*63+1);
            set(pl,'MarkerSize',4,'LineWidth',2,'Color',[ g(c,:) ] )
        end
        colorbar
        colormap(jet)
    end

    if exist('maex', 'var')
        hold on
        pl = plot(maex,-maey,'*m');
        set(pl,'MarkerSize',8,'LineWidth',2)
    end

    if exist('maex', 'var')
        hold on
        pl = plot(maex,-maey,'*m');
        set(pl,'MarkerSize',8,'LineWidth',2)
    end

    xlabel('Distance [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
    ylabel('Depth [km]','FontWeight','bold','FontSize',ZmapGlobal.Data.fontsz.m)
    strib = [  'Misfit '];
    title(strib,'FontWeight','bold',...
        'FontSize',ZmapGlobal.Data.fontsz.m,'Color','k')

    set(gca,'Color',color_bg);
    set(gca,'box','on',...
        'SortMethod','childorder','TickDir','out','FontWeight',...
        'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)
    uicontrol(...
        'Style','pushbutton',...
        'Units','normalized',...
        'Position',[0.9 0.7 0.08 0.08],...
        'String','Grid',...
        'callback',@(~,~)mificrgr(mi,inde));

    uicontrol(...
        'Style','pushbutton',...
        'Units','normalized',...
        'Position',[0.9 0.6 0.08 0.08],...
        'String','Sel EQ',...
        'callback',@cb_pickinv);

    watchoff

    function cb_pickinv(~,~)
        newa2=crossel(newa);
        ZG.newt2=newa2;
        ZG.newcat=newa2;
        timeplot();
    end

end

author: CelsoReyes


text:
function [Index,MAX] = CpuGroup(numberOfGroups,xPrime,numberOfVariables)

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    varsPerGroup = floor(numberOfVariables/numberOfGroups);
    if varsPerGroup == 1
        Index = linspace(1,numberOfVariables,numberOfVariables);
        MAX   = numberOfVariables;
    else
        B      = ones(1,varsPerGroup*numberOfGroups);
        remain = ones(1,(numberOfVariables-varsPerGroup*numberOfGroups))*(numberOfGroups+1);
        R      = reshape(B,varsPerGroup,numberOfGroups);
        k      = linspace(1,numberOfGroups,numberOfGroups);
        index  = R.*repmat(k,varsPerGroup,1);
        index  = reshape(index,1,varsPerGroup*numberOfGroups);
        INDEX  = [index remain];
        [~,I]  = sort(xPrime);
        Index(I) = INDEX;
        if(mod(numberOfVariables,numberOfGroups)==0)
            MAX = numberOfGroups;
        else
            MAX = numberOfGroups + 1;
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function M = volume_intersect_sphere_analytical(varargin)
% Compute the overlap volume between 2 spheres defined in an array.
% Computation is vectorized, and intersection volume are computed an
% analytical way.
%
% Input: spheres data presented in an array G of four columns.
%        G contains parameters of the n spheres
%           . G(1:n,1) - x-coordinate of the center of spheres,
%           . G(1:n,2) - y-coordinate of the center of spheres,
%           . G(1:n,3) - z-coordinate of the center of spheres,
%           . G(1:n,4) - radii of the spheres
%        Each row of the array contains the information for one sphere.
%
%        Input can also be provided in three different vectors. These
%        vectors can be row or column vectors. The 1st one corresponds to
%        x-coordinate of the center of spheres, the 2nd one to the
%        y-coordinate, the 3rd one to the z-coordinate and the 4th one to
%        the radii of the spheres.
%        An optional binary argument can be provided to display or not the
%        result.
%
% Output: Square matrix M(n,n) containing intersection volumes between
% spheres
%         M(i,j) contains the intersection volume between spheres i & j
%         By definition, M(i,i) corresponds to the volume of sphere i.
%
%
% Examples
%
%         x = [0,1,5,3,-5];
%         y = [0,4,3,7,0];
%         z = [0,4,3,7,0];
%         r = [1,5,3,2,2];
%         Display_solution = true;
%         disp('First call')
%         M = volume_intersect_sphere_analytical(x,y,z,r,Display_solution);
%         disp('Second call')
%         M = volume_intersect_sphere_analytical([x',y',z',r'],false);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For 2 spheres i & j, three cases are possible depending on the distance
% d(i,j) of the centers of the spheres i and j.
%   Case 1: spheres i & j do not overlap, there is no overlap volume
%           M(i,j)=0;
%             Condition: d(i,j)>= ri+rj
%             M(i,j) = 0;
%   Case 2: spheres i & j fully overlap, the overlap volume has to be
%           computed.
%             Condition: d(i,j)<= abs(ri-rj)
%             M(i,j) = 4/3*pi*min(ri,rj).^3
%   Case 3: spheres i & j partially overlap, the overlap volume has to be
%           computed decomposing the overlap volume.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Guillaume JACQUENOT
% guillaume dot jacquenot at gmail dot com
% 2008 01 31
% 2009 09 10
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin==0
    % Performs an example
    % Creation of 5 spheres
    x = [0,1,5,3,-5];
    y = [0,4,3,7,0];
    z = [0,4,3,7,0];
    r = [1,5,3,2,2];
    Display_solution = true;
elseif nargin==1 || nargin==2
    temp = varargin{1};
    x = temp(:,1);
    y = temp(:,2);
    z = temp(:,3);
    r = temp(:,4);
    if nargin == 2
        Display_solution = varargin{2};
    else
        Display_solution = false;
    end
elseif nargin==4 || nargin==5
    x = varargin{1};
    y = varargin{2};
    z = varargin{3};
    r = varargin{4};
    if nargin == 5
        Display_solution = varargin{5};
    else
        Display_solution = false;
    end
else
    help volume_intersect_sphere_analytical
    error('volume_intersect_sphere_analytical:e0',...
          'The number of arguments must 0, 1, 2, 4 or 5');
end

% Checking input argument
if ~islogical(Display_solution)
    error('volume_intersect_sphere_analytical:e1',...
          'Display_solution should be logical variable')
end

% Inputs are reshaped
size_x = numel(x);
size_y = numel(y);
size_z = numel(z);
size_r = numel(r);

x = reshape(x,size_x,1);
y = reshape(y,size_y,1);
z = reshape(z,size_z,1);
r = reshape(r,size_r,1);

% Checking if the three input vectors have the same length
if (size_x~=size_y)||(size_x~=size_z)||(size_x~=size_r)
    error('volume_intersect_sphere_analytical:e2',...
          'Input of function must be the same length');
end

% Checking if there is any negative or null radius
if any(r<=0)
    disp(['spheres with null or negative radius'...
          ' won''t be taken into account in the computation.'])
    temp = (r>0);
    x = x(temp);
    y = y(temp);
    z = z(temp);
    r = r(temp);
end

% Checking the size of the input argument
if size_x==1
    M = 4/3*pi*r.^3;
    return
end

% Computation of distance between all spheres, which will allow to
% determine which cases to use.
X = meshgrid(x);
Y = meshgrid(y);
Z = meshgrid(z);
D = sqrt((X-X').^2+(Y-Y').^2+(Z-Z').^2);

% Since the resulting matrix M is symmetric M(i,j)=M(j,i), computations are
% performed only on the upper part of the matrix
D = triu(D,1);

[R1,R2] = meshgrid(r);
sumR    = triu(R1+R2,1);
difR    = triu(abs(R1-R2),1);

% Creating the resulting vector
M = zeros(size_x*size_x,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% C2: Case 2: spheres i & j fully overlap
% One of the spheres is inside the other one.
C2    = triu(D<=difR);
M(C2) = 4/3*pi*min(R1(C2),R2(C2)).^3;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Case 3: spheres i & j partially overlap
% Partial intersection between spheres i & j
C3 = (D>difR)&(D<sumR);
% Computation of the coordinates of one of the intersection points of the
% spheres i & j
Xi(C3,1) = (R1(C3).^2-R2(C3).^2+D(C3).^2)./(2*D(C3));

H1(C3,1) = R1(C3)-Xi(C3,1);
H2(C3,1) = R2(C3)+Xi(C3,1)-D(C3);

% Computation of the partial intersection volume between spheres i & j
M(C3,1) = pi/3*(H1(C3,1).^2.*(3*R1(C3)-H1(C3,1))+...
                H2(C3,1).^2.*(3*R2(C3)-H2(C3,1)));

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

% Compute the volume of each sphere.
% Assign the results to the diagonal of M
M(1:size_x+1:size_x*size_x) = 4/3*pi.*r.^3;

% Conversion of vector M to matrix M
M = reshape(M,size_x,size_x);

% Creating the lower part of the matrix
M = M + tril(M',-1);

% Display results if no argument is provided
if Display_solution
    f = figure;
    hAxs = axes('Parent',f);
    hold on, box on, axis equal
    xlabel('x')
    ylabel('y','Rotation',0)
    zlabel('z','Rotation',0)
    title(['Compute the intersection volume between 2 '...
           ' spheres a vectorized way for several spheres'])
    colour = rand(size_x,3);
    [X,Y,Z]=sphere(20);
    for t = 1: size_x
        surf(x(t)+r(t)*X,y(t)+r(t)*Y,z(t)+r(t)*Z,...
            'FaceColor',colour(t,:),'EdgeColor','none','parent',hAxs)

%         surf(x(t)+r(t)*X,y(t)+r(t)*Y,z(t)+r(t)*Z,...
%             'FaceColor','none','parent',hAxs)
        for u = t+1:size_x
            if M(t,u)~=0
              plot3([x(t),x(u)],[y(t),y(u)],[z(t),z(u)],'k-','parent',hAxs)
              text((x(t)+x(u))/2,(y(t)+y(u))/2,(z(t)+z(u))/2,...
                    num2str(M(t,u)),...
                    'HorizontalAlignment','center',...
                    'BackgroundColor',[1 1 1],'parent',hAxs)
            end
        end
    end
    axis tight
    view([  0.8765    0.7660   -0.0000   -0.8213
           -0.6140    0.3778    0.8090   -0.2864
           -0.8451    0.5200   -0.5878   10.2792
                 0         0         0    1.0000])
end

author: Sable

repo: mcbench-benchmarks


text:
function out = planesBisector(plane1, plane2)
% PLANESBISECTOR  Bisector plane between two other planes
%
%   BIS = planesBisector(PL1, PL2);
%   Returns the planes that contains the intersection line between PL1 and
%   PL2 and that bisect the dihedral angle of PL1 and PL2.
%   Note that computing the bisector of PL2 and PL1 (in that order) returns
%   the same plane but with opposite orientation.
%
%   Example
%     % Draw two planes together with their bisector
%     pl1 = createPlane([3 4 5], [1 2 3]);
%     pl2 = createPlane([3 4 5], [2 -3 4]);
%     % compute bisector
%     bis = planesBisector(pl1, pl2);
%     % setup display
%     figure; hold on; axis([0 10 0 10 0 10]);
%     set(gcf, 'renderer', 'opengl')
%     view(3);
%     % draw the planes
%     drawPlane3d(pl1, 'g');
%     drawPlane3d(pl2, 'g');
%     drawPlane3d(bis, 'b');
%
%   See also
%   planes3d, dihedralAngle, intersectPlanes
%
%   Author: Ben X. Kang
%   Dept. Orthopaedics & Traumatology
%   Li Ka Shing Faculty of Medicine
%   The University of Hong Kong
%   Pok Fu Lam, Hong Kong
%

% Let the two planes be defined by equations
%
%  a1*x + b1*y + c1*z + d1 = 0
%
% and
%
%  a2*x + b2*y + c2*z + d2 = 0
%
% in which vectors [a1,b1,c1] and [a2,b2,c2] are normalized to be of unit
% length (a^2+b^2+c^2 = 1). Then
%
%  (a1+a2)*x + (b1+b2)*y + (c1+c2)*z + (d1+d2) = 0
%
% is the equation of the desired plane which bisects the dihedral angle
% between the two planes.  These coefficients cannot be all zero because
% the two given planes are not parallel.
%
% Notice that there is a second solution to this problem
%
%  (a1-a2)*x + (b1-b2)*y + (c1-c2)*z + (d1-d2) = 0
%
% which also is a valid plane and orthogonal to the first solution. One of
% these planes bisects the acute dihedral angle, and the other the
% supplementary obtuse dihedral angle, between the two given planes.

P1 = plane1(1:3);			% a point on the plane
n1 = planeNormal(plane1);	% the normal of the plane
% d1 = -dot(n1, P1);		% for line equation

P2 = plane2(1:3);
n2 = planeNormal(plane2);
% d2 = -dot(n2, P2);

if ~isequal(P1(1:3), P2(1:3))
	L = intersectPlanes(plane1, plane2);	% intersection of the given two planes
	Pt = L(1:3);							% a point on the line intersection
% 	v2 = cross(n1-n2, L(4:6));				% another vector lie on the bisect plane
% 	out = [v1, v2]';
else
	Pt = P1(1:3);
end

% use column-wise vector
out = createPlane(Pt, n1 - n2);

%%  EOF  %%

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function imsegs = APPimages2superpixels(input_dir, ext, gt)
% imsegs = APPimages2superpixels(input_dir, ext, gt)
% Create the imsegs structure from a segmentation image
%
% INPUT:
% input_dir - The directory of segmentation images. Segments are denoted by
% different colors.
% ext - the extension of the segmentation image filenames
% gt - an existing imsegs structure (or [] if none exists)
% OUTPUT:
% imsegs - image segmentation data (and sometimes ground truth)
%
% Example: imsegs = APPimages2superpixels('../images', 'pnm', [])
%
% Copyright(C) Derek Hoiem, Carnegie Mellon University, 2005
% Permission granted to non-commercial enterprises for
% modification/redistribution under GNU GPL.
% Current Version: 1.0  09/30/2005

if isempty(gt)
    files = dir([input_dir '/*.' ext]);
    for f = 1:length(files)
        gt(f).image_name = files(f).name;
    end
end

imsegs(length(gt)) = struct('imname', '', 'imsize', [0 0]);
for f = 1:length(gt)
    imname = gt(f).image_name;
    disp(imname)
    basename = strtok(imname, '.');
    im = imread([input_dir '/' basename '.' ext]);

    im = double(im);

    imsegs(f).imname = [strtok(imname,'.') '.jpg'];
    imsegs(f).imsize = size(im);
    imsegs(f).imsize = imsegs(f).imsize(1:2);
    im = im(:, :, 1) + im(:, :, 2)*256 + im(:, :, 3)*256^2;
    [gid, gn] = grp2idx(im);
    imsegs(f).segimage = uint16(reshape(gid, imsegs(f).imsize));
    imsegs(f).nseg = length(gn);

end
imsegs = APPgetSpStats(imsegs);

author: zouchuhang

repo: LayoutNet


text:
%% Testing other functions
% Create a Mesh FEM results
clear; close all;

file = 'test2d_triangle';
a.fileName = file;
s = FemDataContainer(a);
mesh = s.mesh;

%% Create functions
% AnalyticalFunction

sAF.fHandle = @(x) x(1,:,:);
sAF.ndimf   = 1;
sAF.mesh    = mesh;
xFun = AnalyticalFunction(sAF);

p1 = xFun.project('P1');
p1.plot

q = Quadrature.create(mesh,'LINEAR');
p1.evaluate(q.posgp)

% get dofs from logical analytical function?
leftCond = @(x) x(:,1) == 0;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [p, stats] = symamd2 (S, knobs)
%SYMAMD Symmetric approximate minimum degree permutation.
%    P = SYMAMD2(S) for a symmetric positive definite matrix S, returns the
%    permutation vector p such that S(p,p) tends to have a sparser Cholesky
%    factor than S.  Sometimes SYMAMD works well for symmetric indefinite
%    matrices too.  The matrix S is assumed to be symmetric; only the
%    strictly lower triangular part is referenced.   S must be square.
%    Note that p = amd(S) is much faster and generates comparable orderings.
%    The ordering is followed by an elimination tree post-ordering.
%
%    Note that this function is source code for the built-in MATLAB symamd
%    function.  It has been renamed here to symamd2 to avoid a filename clash.
%    symamd and symamd2 are identical.
%
%    See also SYMAMD, AMD, COLAMD, COLAMD2.
%
%    Example:
%            P = symamd2 (S)
%            [P, stats] = symamd2 (S, knobs)
%
%    knobs is an optional one- to two-element input vector.  If S is n-by-n,
%    then rows and columns with more than max(16,knobs(1)*sqrt(n)) entries are
%    removed prior to ordering, and ordered last in the output permutation P.
%    No rows/columns are removed if knobs(1)<0.  If knobs(2) is nonzero, stats
%    and knobs are printed.  The default is knobs = [10 0].  Note that knobs
%    differs from earlier versions of symamd.

%    Copyright 1998-2007, Timothy A. Davis, and Stefan Larimore
%    Developed in collaboration with J. Gilbert and E. Ng.
%    Acknowledgements: This work was supported by the National Science
%       Foundation, under grants DMS-9504974 and DMS-9803599.

%-------------------------------------------------------------------------------
% perform the symamd ordering:
%-------------------------------------------------------------------------------

if (nargout <= 1 & nargin == 1)						    %#ok
    p = symamd2mex (S) ;
elseif (nargout <= 1 & nargin == 2)					    %#ok
    p = symamd2mex (S, knobs) ;
elseif (nargout == 2 & nargin == 1)					    %#ok
    [p, stats] = symamd2mex (S) ;
elseif (nargout == 2 & nargin == 2)					    %#ok
    [p, stats] = symamd2mex (S, knobs) ;
else
    error('symamd:  incorrect number of input and/or output arguments.') ;
end

%-------------------------------------------------------------------------------
% symmetric elimination tree post-ordering:
%-------------------------------------------------------------------------------

[ignore, q] = etree (S (p,p)) ;
p = p (q) ;

%    stats is an optional 20-element output vector that provides data about the
%    ordering and the validity of the input matrix S.  Ordering statistics are
%    in stats (1:3).  stats (1) = stats (2) is the number of dense or empty
%    rows and columns ignored by SYMAMD and stats (3) is the number of
%    garbage collections performed on the internal data structure used by
%    SYMAMD (roughly of size 8.4*nnz(tril(S,-1)) + 9*n integers).
%
%    MATLAB built-in functions are intended to generate valid sparse matrices,
%    with no duplicate entries, with ascending row indices of the nonzeros
%    in each column, with a non-negative number of entries in each column (!)
%    and so on.  If a matrix is invalid, then SYMAMD may or may not be able
%    to continue.  If there are duplicate entries (a row index appears two or
%    more times in the same column) or if the row indices in a column are out
%    of order, then SYMAMD can correct these errors by ignoring the duplicate
%    entries and sorting each column of its internal copy of the matrix S (the
%    input matrix S is not repaired, however).  If a matrix is invalid in other
%    ways then SYMAMD cannot continue, an error message is printed, and no
%    output arguments (P or stats) are returned.  SYMAMD is thus a simple way
%    to check a sparse matrix to see if it's valid.
%
%    stats (4:7) provide information if SYMAMD was able to continue.  The
%    matrix is OK if stats (4) is zero, or 1 if invalid.  stats (5) is the
%    rightmost column index that is unsorted or contains duplicate entries,
%    or zero if no such column exists.  stats (6) is the last seen duplicate
%    or out-of-order row index in the column index given by stats (5), or zero
%    if no such row index exists.  stats (7) is the number of duplicate or
%    out-of-order row indices.
%
%    stats (8:20) is always zero in the current version of SYMAMD (reserved
%    for future use).

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% test the spherical compression scheme

sampling_type = 'random';
% sampling_type = 'uniform';

data_type = 'syntetic';
% data_type = 'real';

n = 1500;
k = 2;          % number of vanishing moments.
dir = 1;

switch sampling_type
    case 'random'
        pos = randn(3,n);
        d = sqrt( pos(1,:).^2+pos(2,:).^2+pos(3,:).^2 );
        pos(1,:) = pos(1,:)./d;
        pos(2,:) = pos(2,:)./d;
        pos(3,:) = pos(3,:)./d;
    case 'uniform'
        [pos,L,diam] =  sphere_sampling(n);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute some function
switch data_type
    case 'syntetic'
        v = pos(1,:).^2 - pos(2,:).*pos(3,:);
        v = v(:);
        v = rescale(v);
        v = cos(10*v);
        % list for compression plot
        m_list = [1,2,5,10,15];
    case 'real'
        sampling_type = 'real';
        load('data/data_lenglet.mat');
%        [pos,face] = gen_base_mesh('oct');
%        [pos,face] = subdivide_sphere(pos,face,2); pos = pos';
        v = ADC24_36_14(:);
        pos = grad_81';
        % symmetrize the data
        pos = [pos,-pos];
        v = [v;v];
        n = length(v);
        v = rescale(v);
        % list for compression plot
        m_list = [5,10,20,40,70];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute clustering
clear options;
options.ptmax = k^2;
options.part_type = 'kmeans';
[part,B,G] = dichotomic_grouping(pos,options);
% plot
clf;
subplot(1,2,1);
plot_spherical_partition(pos,part);
title('Clustering of sampling locations');

subplot(1,2,2);
plot_spherical_function(pos,v);
title('Function to transform (interpolated)');
saveas(gcf, ['function_',data_type, '_clustering_', sampling_type], 'png');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% perform transform
clear options;
options.part_type = 'kmeans';
[w,info,part] = perform_alpert_transform_nd(v,pos,k, 1, options);
% try reconstruction just to check the transform is bijective ...
options.part = part;
v1 = perform_alpert_transform_nd(w,pos,k, -1,options);
norm(v-v1, 'fro') % should be 0

cv = rev_sort_abs(v);
cw = rev_sort_abs(w);
cv = l2error(v);
cw = l2error(w);

a = (1:n/2)';
x = log2(a);
y1 = cv(a);
y2 = cw(a);

% reglin
fiton = a;
A = a;
B = y2;
aa = log2(A(fiton));
bb = log2(B(fiton));
[coef,S] = polyfit(aa,bb,1);
disp( sprintf('----> 1D Alpert transform, decreasing=%.4f', coef(1)) );
reglin = log2(a)*coef(1)+coef(2);

clf;
hold on;
plot( x, log2(y1), x, log2(y2), x, reglin, 'b:' );
axis tight;
legend('Original error decreasing', 'Transformed error decreasing');
hold off;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% perform reconstruction
clf;
subplot(2,3,1);
plot_spherical_function(pos,v);
title('Original function');
s = 1;

for m = m_list
    s = s+1;
    w1 = keep_biggest(w, round(n*m/100) );
    % backward transform
    v1 = perform_alpert_transform_nd(w1,pos,k, -1,part);
    subplot(2,3,s);
    plot_spherical_function(pos,v1);
    err = round( 100*norm(v-v1, 'fro')^2/norm(v, 'fro')^2 );    % RMS error in %
    title(sprintf('%.1f %% of coef, err=%d%%', m, err  ));
end

saveas(gcf, ['compression_',data_type, '_clustering_', sampling_type], 'png');

author: gpeyre

repo: matlab-toolboxes


text:
%RANDRESET  Reset state of random generators
%
%   RANDRESET                   - reset states to 1
%   RANDRESET(STATE)            - reset states to STATE
%   STATE = RANDRESET           - retrieve present state (no reset)
%   OLDSTATE = RANDRESET(STATE) - retrieve present state and reset

function output = randreset(state)

if nargin < 1 & nargout < 1
   state = 1;
end

if nargin > 0
  if iscell(state)
    % old state format, force switch to old generator
    rand('seed',0);
    % new state format, force switch to new generator
  elseif isstruct(state)
    rng(0,'twister');
  end
end

if verLessThan('matlab','8.1.0') || strcmp(getfield(rng,'Type'),'Legacy')
  % old random number generators

  if nargout > 0
    % retrieve present state
    randstate = cell(1,2);
    randstate{1} = rand('state');
    randstate{2} = randn('state');
    output = randstate;
  end

  if ~(nargout == 1 & nargin == 0)
    % reset
    if iscell(state)
      rand('state',state{1});
      randn('state',state{2});
    else
      rand('state',state);
      randn('state',state);
    end
  end

else
  % new random number generator (rng)

  if nargout > 0
    % retrieve present state
    output = rng;
  end

  if ~(nargout == 1 & nargin == 0)
    % reset
    rng(state);
  end

end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function varargout = setStateName(obj,varargin)
    nVarargs = length(varargin);
    name = varargin{1};
    if nVarargs == 2
        index = varargin{2};
    else
        index = 1:obj.dim.x;
    end
    varargout = cell(size(index));
    j = 1;
    for i=index
        obj.x(i) = sym(name{j});
        varargout{j} = sym(name{j});
        j = j+1;
    end
end

author: deng-haoyang

repo: ParNMPC


text:
function [R,lambda] = linreg( X , Y , lambda )

%X = [ones(size(X,1),1) X];

%% method 1: soving linear regression using close-form solution %%%%%%%%%%%

%R = (X'*X+eye(size(X,2))*lambda)\X'*Y;

%% method 2: using SVR in liblinear %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
featdim   = size(X,2);
shapedim  = size(Y,2);

param = sprintf('-s 12 -p 0 -c %f -q', lambda);
%param = sprintf('-s 12 -p 0 -c 0.3 -q');
R_tmp = zeros( featdim, shapedim );
tic;
for o = 1 : shapedim
    disp(['Training landmarks ' num2str(o)]);
    model = train(Y(:,o),sparse(X),param);
    R_tmp(:,o) = model.w';
end
toc;

R = R_tmp;

end

author: tntrung


text:
function Align = Find_Fiducials(q, r, s, Align, ecg, Param)
%OVERVIEW, Finds the lead with maximum QRS amplitude, takes a beat as QRS template and finds fiducial points such that
% aligned beats have maximum cross-correlation (over Q to S or S to T interval) in the chosen
% lead against the template.
%
% INPUTS            MANDATORY           OUTPUTS
%                   q                   q point locations
%
%                   r                   r point locations
%
%                   s                   s point locations
%
%                   Align               Alignment information for beats in
%                                       analysis window.
%
%                   ecg                 N by M array of ecg data. Contains M channels of ecg with N data points each.
%
%                   Param               Structure containing the parameters
%                                       used for evaluating TWAs.
%
% OUTPUTS
%                   Align               Alignment information used to
%                                       determine if current window is clean enough for
%                                       analysis.
%
%   REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Written by Shamim Nemati,
%       editted by Ismail Sadiq on 10/26/2019.
%	COPYRIGHT (C) 2019
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information. The license may be found in
%       the Documents folder of the Physionet-Cardiovascular-Signal-Toolbox.

q = q(1:min(length(q), length(s)));
s = s(1:min(length(q), length(s)));

% Sort leads by amplitude.
leadSorted = Sort_Leads(q, s, Align.st, ecg, Param);

for ilead = 1:length(leadSorted)
    Align.lead = leadSorted(ilead);

    % Get the orientation of the QRS complex, max positive or max negative deflection.
    Align.orientation = MaxFrontOrientation(ecg(:, Align.lead), q, s);

    Align.fidBase = r;

    Align.f2s = floor(median(s - Align.fidBase) + Param.stAdjIntv);   % reserve to avoid including part of QRS in ST interval computations
    Align.q2f = min([floor(median(Align.fidBase - q)),Align.f2s-Param.stAdjIntv]);

    Align.fid = [];
    NFailed = 0;

    Align.template = 0;
    while isempty(Align.fid)
        Align.template = Align.template + 1;
        % Get alignment information
        Align = Adjust_e(Align, ecg, 1000, Param);
        if isempty(Align.fid)

            NFailed = NFailed + 1;
            if (NFailed > 0.1 * length(s))
                break
            end
        end
    end
    if (~isempty(Align.fid))
        for lead=1:size(ecg, 2)
            for i = 1:length(Align.fid)
                Align.amp(i, lead) = (max(ecg(q(i) : s(i), lead))+min(ecg(q(i) : s(i), lead)))/2;  % align to zero mean on QS segment
            end
        end
        % get correlation for each beat with template over Q-S and S-T
        % interval
        Align = CheckLead_Correlations(Align, ecg, Param);
        disp(['FindFiducials: Alignment in lead ' num2str(Align.lead)]);
        break;
    end;
    disp(['FindFiducials: lead ' num2str(Align.lead) ' is too noisy or ectopic..']);
end;
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function leads = Sort_Leads(q, s, st, ecg, Param)

NLead = size(ecg, 2);

%   QRS or T amplitudes in each lead
ptp = zeros(min(length(q), 8),NLead);
for i = 1:min(length(q), 8)
    for l = 1:NLead
        if strcmp(Param.Alignment, 'st')
            ptp(i, l) = max(ecg(s(i):s(i) + st, l)) - min(ecg(s(i):s(i) + st, l));
        else
            ptp(i, l) = max(ecg(q(i):s(i), l)) - min(ecg(q(i):s(i), l));
        end;
    end;
end;
ptpmed = zeros(1,NLead);
for i = 1:NLead
    ptpmed(i) = median(ptp(:, i));
end;

[sorted indices] = sort(ptpmed);
leads = indices(length(indices):-1:1);
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Align = Adjust_e(Align, ecg, freq, Param)
% OVERVIEW, This function determines which beats in analysis window for
% the lead given by Align.lead are valid for analysis, based on correlation
% with a QRS and S-T offset segment template.
Align.fidBase(Align.template);

len_ecg = size(ecg,1);
qs_templ = ecg(max(Align.fidBase(Align.template) - Align.q2f,1):min(Align.fidBase(Align.template) + Align.f2s,len_ecg), Align.lead);
st_templ = ecg(max(Align.fidBase(Align.template) + Align.f2s,1):min(Align.fidBase(Align.template) + Align.f2s + Align.st,len_ecg), Align.lead);

[a b c d] = Adjust_Fiducials(ecg, Align.fidBase, Align.q2f, Align.f2s, qs_templ, Align.st, st_templ, Param);
% fid a -> with respect to QRS alignment
% fid c -> with respect to ST alignment

if ~isempty(c)
    Align.fidQRS = a;
    Align.QRScorr(:, Align.lead) = b;
    Align.fid = c;
    Align.Tcorr(:, Align.lead) = d;
    amp = (max(qs_templ)+min(qs_templ))/2;
    Align.qs_templ(:, Align.lead) = qs_templ - amp;
    Align.st_templ(:, Align.lead) = st_templ - amp;

    for i = 1:length(Align.fid)
        Align.valid(i, Align.lead) = (Align.QRScorr(i, Align.lead) >= Param.corrQRS) && (Align.Tcorr(i, Align.lead) >= Param.corrT);
    end;
    Align.validleads(Align.lead) = true;
end;

return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function orientation = MaxFrontOrientation(ecg, q, s)

NPos = 0;

for i = 1:min(20, length(s))
    [ma maI] = max(ecg(q(i):s(i)));
    [mi miI] = min(ecg(q(i):s(i)));

    NPos = NPos + (maI > miI);
end;

orientation = (NPos > 10);

return;

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function J = computeCostMulti(X, y, theta)
%COMPUTECOSTMULTI Compute cost for linear regression with multiple variables
%   J = COMPUTECOSTMULTI(X, y, theta) computes the cost of using theta as the
%   parameter for linear regression to fit the data points in X and y

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta
%               You should set J to the cost.

J = (X*theta - y)' * (X*theta - y) / (2*m);

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

end

author: merwan

repo: ml-class

save_path: github-repos/MATLAB/merwan-ml-class


text:
function [detjacobian,invjacobian]=Jacobian(nnel,dhdr,dhds,xcoord,ycoord)

%------------------------------------------------------------------------
%  Purpose:
%     determine the Jacobian for two-dimensional mapping
%
%  Synopsis:
%     [detjacobian,invjacobian]=Jacobian(nnel,dhdr,dhds,xcoord,ycoord)
%
%  Variable Description:
%     jacobian - Jacobian for one-dimension
%     nnel - number of nodes per element
%     dhdr - derivative of shape functions w.r.t. natural coordinate r
%     dhds - derivative of shape functions w.r.t. natural coordinate s
%     xcoord - x axis coordinate values of nodes
%     ycoord - y axis coordinate values of nodes
%------------------------------------------------------------------------

 jacobian=zeros(2,2);

 for i=1:nnel
 jacobian(1,1) = jacobian(1,1)+dhdr(i)*xcoord(i);
 jacobian(1,2) = jacobian(1,2)+dhdr(i)*ycoord(i);
 jacobian(2,1) = jacobian(2,1)+dhds(i)*xcoord(i);
 jacobian(2,2) = jacobian(2,2)+dhds(i)*ycoord(i);
 end

 detjacobian = det(jacobian) ;  % Determinant of Jacobian matrix
 invjacobian = inv(jacobian) ;  % Inverse of Jacobian matrix


text:
function res= experiment_run(exp,net)

net.cnn_id = exp.model;
switch net.cnn_mode
case 0 % matconvnet
    exp.opts.normalize = get_cnn_normalize(net.normalization) ;
    exp.opts.denormalize = get_cnn_denormalize(net.normalization) ;
case 1
    switch exp.model
    case {'alex','nin','gnet'} % alexnet, nin, inception
        d = load('data/ilsvrc_2012_mean.mat');
        mm = d.image_mean;
    case 'vgg16' % vgg
        mm = repmat(reshape([103.939, 116.779, 123.68],[1,1,3]),[256 256]);
    case 'alexWeb'
        mm = repmat(reshape([0.8309, 0.8310, 0.8328],[1,1,3]),[256 256]);
    end
    net.normalization=struct('imageSize',[net.im_sz net.im_sz 3],'averageImage',mm);
    exp.opts.normalize = @(x) U_prepare_image(x, mm,exp.model,-2);
    exp.opts.denormalize = @(x) U_prepare_image(x, mm,exp.model,-1);
end

net.task = exp.opts.task;
switch exp.opts.task
case {1,2}
    % 1: neuron inversion
    % 2: neuron inpainting
    res = invert_nn_dw(net, exp.opts.mask, exp.opts) ;
case 0
    % feature inversion
    % same as deep-goggle
    if ischar(exp.opts.feats)
        im = single(imread(exp.opts.feats));
        if size(im,3) == 1, im = cat(3,im,im,im) ; end
        exp.opts.feats = compute_features(net, exp.opts.normalize(im),0);
        %keyboard
    end
    % gradient-descent
    net.feats = exp.opts.feats;
    res = invert_nn_dw(net, exp.opts.feats, exp.opts) ;
end

end
function args = expandOpts(opts)
% -------------------------------------------------------------------------
args = horzcat(fieldnames(opts), struct2cell(opts))' ;
end

author: donglaiw

repo: mNeuron

save_path: github-repos/MATLAB/donglaiw-mNeuron


text:
%% file example_mtl_classify.m
%   This example shows how to perform classification using least squares
%   loss and logistic loss.
%
%% LICENSE
%   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/>.
%
%   Copyright (C) 2011 - 2012 Jiayu Zhou and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 3, 2012.
%

clear, clc;
addpath('../MALSAR/utils/')
addpath('../MALSAR/functions/Lasso/')

n = 50;
d = 300;
t = 10;

X = cell(t, 1);
Y = cell(t, 1);
W = randn(d, t);
W_mask = abs(randn(d, t))<1;
W(W_mask) = 0;
for i = 1: t
    X{i} = randn(n, d);
    Y{i} = sign(X{i} * W(:, i) + rand(n, 1) * 0.01);
end

% training and prediction using least squares loss
W_pred = Least_Lasso(X, Y, 0.01);
% compute training error
least_acc = zeros(t, 1);
for i = 1: t
    least_acc(i) = nnz(sign(X{i} * W_pred(:, i)) == Y{i})/n;
end
fprintf('Least Squares Loss Training Accuracy: %.4f +/- %.4f\n', mean(least_acc), std(least_acc));

% training and prediction using logistic loss
[W_pred C_pred]= Logistic_Lasso(X, Y, 0.01);
% compute training error
logistic_acc = zeros(t, 1);
for i = 1: t
    logistic_acc(i) = nnz(sign(X{i} * W_pred(:, i) + C_pred(i)) == Y{i})/n;
end
fprintf('Logistic Loss Training Accuracy: %.4f +/- %.4f\n', mean(logistic_acc), std(logistic_acc));

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function variance = chi_square_noncentral_variance ( a, b )

%*****************************************************************************80
%
%% CHI_SQUARE_NONCENTRAL_VARIANCE returns the variance of the noncentral Chi squared PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, the parameter of the PDF.
%    1 <= A.
%
%    Input, real B, the noncentrality parameter of the PDF.
%    0.0 <= B.
%
%    Output, real VARIANCE, the variance value.
%
  variance = 2.0 * ( a + 2.0 * b );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function DEM_demo_DCM_LAP
% Demo applying the Laplace scheme to DCM with hidden states
%__________________________________________________________________________
% This routine demonstrates Generalized filtering for a DCM (Dynamic Causal
% Model) of fMRI responses using simulated data. This is an endogenous
% DCM in that there are no exogenous inputs. The demonstration specifies
% and inverts a full connectivity model and then illustrates post-hoc model
% optimization to recover (discover) the true architecture. It concludes
% with an automatic model optimization in terms of the prior variances over
% coupling parameters.
%__________________________________________________________________________
% Copyright (C) 2010 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: DEM_demo_DCM_LAP.m 6483 2015-06-21 21:14:34Z karl $

% Specify a DCM to generate synthetic data
%==========================================================================
rng('default')

% DEM Structure: create inputs
% -------------------------------------------------------------------------
T  = 256;
TR = 3.22;
n  = 4;
U  = spm_conv(randn(n,T),0,2)/4;

% set inversion parameters
% -------------------------------------------------------------------------
M(1).E.s  = 1/2;         % smoothness of random fluctuations
M(1).E.d  = 2;           % embedding dimension
M(1).E.n  = 4;           % embedding dimension
M(1).E.nE = 32;          % maximum number of DEM iterations

% priors
% -------------------------------------------------------------------------
A  = ones(n,n);
B  = zeros(n,n,0);
C  = zeros(n,n);
D  = zeros(n,n,0);

[pE,pC] = spm_dcm_fmri_priors(A,B,C,D);

% adjust M.f (GF works in time bins not seconds) and initialize M.P
% -------------------------------------------------------------------------
M(1).f  = inline(['spm_fx_fmri(x,v,P)*' num2str(TR)],'x','v','P');
M(1).g  = 'spm_gx_fmri';
M(1).x  = sparse(n,5);
M(1).pE = pE;
M(1).pC = pC;

M(2).v  = sparse(n,1);

% simulate endogenous dynamics
%==========================================================================

% true parameters (stochastic estimates)
% -------------------------------------------------------------------------
pP   = pE;
a    =  0.3;
b    = -0.3;
c    =  0.0;
pP.A = [c  a  0  0  0;
        a  c  b  0  0;
        0  b  c  a  0;
        0  0  a  c  b;
        0  0  0  a  c];

pP.A = pP.A(1:n,1:n);
pP.C = eye(n,n);
SIM  = spm_DEM_generate(M,U,{pP},{6,16},{16});

% Show simulated response
%--------------------------------------------------------------------------
spm_figure('Getwin','Figure 1');
spm_DEM_qU(SIM.pU)

% Specify generative model for inversion (DCM)
% =========================================================================

% set inversion parameters
% -------------------------------------------------------------------------
DCM.M       = M;
DCM.M(2).v  = 0;

% allow (only) neuronal [x, s, f, q, v] hidden states to fluctuate
% -------------------------------------------------------------------------
W           = ones(n,1)*exp([12 16 16 16 16]);
DCM.M(1).xP = exp(6);
DCM.M(1).V  = exp(6);        % prior log precision (noise)
DCM.M(1).W  = diag(W);       % fixed precision (hidden-state)
DCM.M(2).V  = exp(16);       % fixed precision (hidden-cause)

% Add data
% -------------------------------------------------------------------------
DCM.Y       = SIM.Y;

% Full connectivity inversion
% =========================================================================
F  = ones(n,n);
B  = zeros(n,n,0);
C  = zeros(n,1);
D  = zeros(n,n,0);

[pE pC]     = spm_dcm_fmri_priors(F,B,C,D);
DCM.M(1).pE = pE;
DCM.M(1).pC = pC;
FULL        = spm_LAP(DCM);

% Search model space with Savage-Dickey density ratio
% =========================================================================
[A K Nk]  = spm_dcm_sparse_priors(n);

% find true model
% -------------------------------------------------------------------------
for i = 1:length(A);
    if ~any(spm_vec(~~(pP.A + eye(n,n)) - A{i}))
        tA = i; break
    end
end

% find candidate models based on full-connectivity
% -------------------------------------------------------------------------
pE    = FULL.M(1).pE.A;
qE    = FULL.qP.P{1}.A;
qC    = FULL.qP.C(1:n*n,1:n*n);
pC    = FULL.M(1).pC(1:n*n,1:n*n);
for i = 1:length(A)
    k       = find(~A{i});
    rE      = pE;
    rC      = pC;
    rE(k)   = 0;
    rC(k,k) = 0;
    P(i,1)  = spm_log_evidence(qE,qC,pE,pC,rE,rC);
end

% posterior density under best model
% -------------------------------------------------------------------------
[p,i]     = max(P);
k         = find(~A{i});
rE        = pE;
rC        = pC;
rE(k)     = 0;
rC(k,k)   = 0;
[F,sE,sC] = spm_log_evidence_reduce(qE,qC,pE,pC,rE,rC);

% log-posterior (model)
% -------------------------------------------------------------------------
PP    = exp(P - max(P));
PP    = PP/sum(PP);

% Graphics (density on parameter and model space)
% -------------------------------------------------------------------------
spm_figure('Getwin','Figure 2'); clf

subplot(2,2,1)
spm_plot_ci(spm_vec(sE),sC),     hold on
bar(spm_vec(pP.A),1/2), hold off
title('true and MAP connections','FontSize',16)
axis square

subplot(2,2,2)
bar(P)
title('log-posterior','FontSize',16)
xlabel('model','FontSize',12)
ylabel('log-probability','FontSize',12)
axis square

subplot(2,2,3)
plot(Nk,    P,    '.k','MarkerSize',16), hold on
plot(Nk(tA),P(tA),'.r','MarkerSize',32), hold off
title('log-evidence','FontSize',16)
xlabel('graph size','FontSize',12)
ylabel('log-probability','FontSize',12)
axis square

subplot(2,2,4)
bar(PP)
title('posterior','FontSize',16)
xlabel('model','FontSize',12)
ylabel('probability','FontSize',12)
axis square

% Compare true and AMS adjacency
% -------------------------------------------------------------------------
spm_figure('Getwin','Figure 3'); clf

% Automatic model selection (optimizing prior variances on parameters)
% =========================================================================
[rE,rC]   = spm_dcm_optimise(qE,qC,pE,pC);
rA        = spm_unvec(diag(rC),pE);

subplot(2,2,1)
imagesc(full(A{tA} - diag(diag(A{tA}))))
title('true adjacency','FontSize',16)
xlabel('source','FontSize',12)
ylabel('target','FontSize',12)
axis square

subplot(2,2,2)
imagesc(full(rA))
title('optmised priors','FontSize',16)
xlabel('source','FontSize',12)
ylabel('target','FontSize',12)
axis square

author: spm


text:
%% Machine Learning Online Class - Exercise 2: Logistic Regression
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the second part
%  of the exercise which covers regularization with logistic regression.
%
%  You will need to complete the following functions in this exericse:
%
%     sigmoid.m
%     costFunction.m
%     predict.m
%     costFunctionReg.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% Load Data
%  The first two columns contains the X values and the third column
%  contains the label (y).

data = load('ex2data2.txt');
X = data(:, [1, 2]); y = data(:, 3);

plotData(X, y);

% Put some labels
hold on;

% Labels and Legend
xlabel('Microchip Test 1')
ylabel('Microchip Test 2')

% Specified in plot order
legend('y = 1', 'y = 0')
hold off;

%% =========== Part 1: Regularized Logistic Regression ============
%  In this part, you are given a dataset with data points that are not
%  linearly separable. However, you would still like to use logistic
%  regression to classify the data points.
%
%  To do so, you introduce more features to use -- in particular, you add
%  polynomial features to our data matrix (similar to polynomial
%  regression).
%

% Add Polynomial Features

% Note that mapFeature also adds a column of ones for us, so the intercept
% term is handled
X = mapFeature(X(:,1), X(:,2));

% Initialize fitting parameters
initial_theta = zeros(size(X, 2), 1);

% Set regularization parameter lambda to 1
lambda = 1;

% Compute and display initial cost and gradient for regularized logistic
% regression
[cost, grad] = costFunctionReg(initial_theta, X, y, lambda);

fprintf('Cost at initial theta (zeros): %f\n', cost);

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ============= Part 2: Regularization and Accuracies =============
%  Optional Exercise:
%  In this part, you will get to try different values of lambda and
%  see how regularization affects the decision coundart
%
%  Try the following values of lambda (0, 1, 10, 100).
%
%  How does the decision boundary change when you vary lambda? How does
%  the training set accuracy vary?
%

% Initialize fitting parameters
initial_theta = zeros(size(X, 2), 1);

% Set regularization parameter lambda to 1 (you should vary this)
lambda = 0;

% Set Options
options = optimset('GradObj', 'on', 'MaxIter', 400);

% Optimize
[theta, J, exit_flag] = ...
	fminunc(@(t)(costFunctionReg(t, X, y, lambda)), initial_theta, options);

% Plot Boundary
plotDecisionBoundary(theta, X, y);
hold on;
title(sprintf('lambda = %g', lambda))

% Labels and Legend
xlabel('Microchip Test 1')
ylabel('Microchip Test 2')

legend('y = 1', 'y = 0', 'Decision boundary')
hold off;

% Compute accuracy on our training set
p = predict(theta, X);

fprintf('Train Accuracy: %f\n', mean(double(p == y)) * 100);

author: Borye


text:
function [test_failed, atomsArr]= test_libltfat_dgtrealmp(varargin)
test_failed = 0;
doplot = 0;

fprintf(' ===============  %s ================ \n',upper(mfilename));

definput.flags.complexity={'double','single'};
[flags]=ltfatarghelper({},definput,varargin);
dataPtr = [flags.complexity, 'Ptr'];
dataPtrPtr = [flags.complexity, 'PtrPtr'];

intbitsize = 8*calllib('libltfat','ltfat_int_size');
intPtr = sprintf('int%dPtr',intbitsize);

[~,~,enuminfo]=libltfatprotofile;
phaseconv = enuminfo.ltfat_phaseconvention;
hintstruct = enuminfo.ltfat_dgt_hint;

algmpstruct = enuminfo.ltfat_dgtmp_alg;
statusenum = enuminfo.ltfat_dgtmp_status;

fftwflags = struct('FFTW_MEASURE',0,'FFTW_ESTIMATE',64,'FFTW_PATIENT',32,'FFTW_DESTROY_INPUT',1,...
    'FFTW_UNALIGNED',2,'FFTW_EXHAUSTIVE',8,'FFTW_PRESERVE_INPUT',16);

base = 2048;
Larr  = [128* base   9   2];
glarr = [ base  10   9   1];
aarr  = [   base/4   10   3   1];
Marr  = [ base  36   3   2];
Warr  = [  1   3   3   1];
errtoldb = -40;

atomsArr= [];

for idx = 1:1%numel(Larr)

%     L = Larr(idx);
%     W = Warr(idx);
%     a = aarr(idx);
%     M = Marr(idx);
%     M2 = floor(M/2) + 1;
%     gl = glarr(idx);
% f = randn(L,W,flags.complexity);

for ii=  57

%filename = sprintf('~/Desktop/SQAM/%02d.wav',ii);
%disp(filename);

%[f,fs] = wavload(filename);
[f,fs] = gspi;

%f = postpad(f,fs);
f = cast(f,flags.complexity);
f = f(:,1);

[Ls,W] = size(f(:,1));
%Ls = min([Ls,2*fs]);
% Ls = 5*2048;
%Ls = 20000;
% f = postpad(f(fs:fs+Ls),Ls);
% f = zeros(Ls,1);
% f(1) = 1;
%f = [zeros(10*1024,1);f;zeros(10*1024,1);];

Ls = numel(f);
%f(:) = pconv(f,fir2long(firwin('hann',1),Ls));
%f = [zeros(100,1);f(6*fs+1:7*fs)];

%f = pconv(f,fir2long(firwin('hann',1),numel(f)));
%f = postpad(f,Ls+fs);
%f = [zeros(fs,1);postpad(f(fs:end,1),Ls);zeros(fs,1)];
W = 1;
%Ls = numel(f);
a  = [     1024,   128 ,   256, 128,  64,  1024,   512, 16, 256];
M  = [     4096,  2048,  1024, 512,  256, 2048, 32, 512];
gl = [     4096,  2048,  1024, 512,  256, 2048,   2048,  32, 512 ];
M2 = floor(M/2) + 1;
P = [2];
Psize = numel(P);

L = dgtlength(Ls,max(a(P)),max(M(P)));
%L = dgtlength(Ls,a(1),M(1));
f = postpad(f,L);

%f(2:2:end) = -f(2:2:end);
cphaseconv = phaseconv.LTFAT_TIMEINV;
mphaseconv = 'timeinv';
if cphaseconv == phaseconv.LTFAT_FREQINV
    mphaseconv = 'freqinv';
end
%f(:) = linspace(0,1,L);

N = L./a;
%g = randn(gl,1,flags.complexity);

gCell = cell(Psize,1);
for p=1:Psize
    gCell{p} = cast(firwin('blackman',gl(P(p)),'peak'),flags.complexity);
%     if p==1
%        gg = cast(firwin('tria',gl(P(p))/2,'peak'),flags.complexity);
%        gCell{p}(end/2 +1:end) = 0;
%        gCell{p}(end- numel(gg(end/2+1:end)) +1:end) = gg(end/2+1:end);
%     end
    gCell{p} = normalize(gCell{p},'2');
%        [A,B] = gabframebounds(gCell{p},a(P(p)),M(P(p)));
%     gCell{p} = cast(gabwin({'gauss',a(p)*M(p)/L},a(p),M(p),L),flags.complexity);
%     [idx] = find(gCell{p} < 1e-6 ,1,'first');
%     gCell{p} = middlepad(gCell{p},2*idx-1);
%     gl(p) = numel(gCell{p});
end
g = cell2mat(gCell);

gPtr = libpointer(dataPtr,g);
glPtr = libpointer(intPtr,gl(P));
aPtr  = libpointer(intPtr,a(P));
MPtr  = libpointer(intPtr,M(P));

fPtr = libpointer(dataPtr,f);

fout = randn(L,numel(P),flags.complexity);
foutPtr = libpointer(dataPtr,fout);

sizeaccum = 0;
for p=1:Psize
    sizeaccum = sizeaccum + M2(P(p))*N(P(p))*W;
end

cout = complex2interleaved(...
cast(zeros(sizeaccum,1)+...
         1i*zeros(sizeaccum,1),flags.complexity));
coutPtr = libpointer(dataPtr,cout);

%ctrue = dgt(f,g(1:gl(1)),a(1),M(1));
atoms = sizeaccum*10;
%atoms = 13;

tic;
params = calllib('libltfat','ltfat_dgtmp_params_allocdef');
calllib('libltfat','ltfat_dgtmp_setpar_maxatoms',params,atoms);
calllib('libltfat','ltfat_dgtmp_setpar_maxit',params,2*atoms);
calllib('libltfat','ltfat_dgtmp_setpar_errtoldb',params,errtoldb);
calllib('libltfat','ltfat_dgtmp_setpar_kernrelthr',params,1e-4);
calllib('libltfat','ltfat_dgtmp_setpar_phaseconv',params,cphaseconv);
calllib('libltfat','ltfat_dgtmp_setpar_pedanticsearch',params,1);
calllib('libltfat','ltfat_dgtmp_setpar_alg',params,algmpstruct.ltfat_dgtmp_alg_loccyclicmp);
%calllib('libltfat','ltfat_dgtmp_setpar_alg',params,algmpstruct.ltfat_dgtmp_alg_LocOMP);
calllib('libltfat','ltfat_dgtmp_setpar_iterstep',params,1e6);
calllib('libltfat','ltfat_dgtmp_setpar_cycles',params,1);

plan = libpointer();
funname = makelibraryname('dgtrealmp_init_gen_compact',flags.complexity,0);
statusInit = calllib('libltfat',funname,gPtr,glPtr,...
    L,Psize,aPtr,MPtr,params,plan);
tinit = toc;

calllib('libltfat','ltfat_dgtmp_params_free',params);

tic
funname = makelibraryname('dgtrealmp_reset',flags.complexity,0);
statusReset = calllib('libltfat',funname,plan,fPtr);
t1 = toc;

 cres1 = complex2interleaved(...
 cast(randn(sizeaccum,1)+...
          1i*randn(sizeaccum,1),flags.complexity));
cresPtr = libpointer(dataPtr,cres1);

if doplot
funname = makelibraryname('dgtrealmp_getresidualcoef_compact',flags.complexity,0);
calllib('libltfat',funname,plan,cresPtr);
cres2 = reshape(postpad(interleaved2complex(cresPtr.value),M2(P(1))*N(P(1))),M2(P(1)),N(P(1)));
figure(2); plotdgtreal(cres2,1,100,'clim',[-90,10]);
end

tic
 funname = makelibraryname('dgtrealmp_execute_compact',flags.complexity,0);
 statusExecute = calllib('libltfat',funname,plan,fPtr,coutPtr,foutPtr);
 %funname = makelibraryname('dgtrealmp_execute_niters_compact',flags.complexity,0);
 %statusExecute = calllib('libltfat',funname,plan,50*atoms,coutPtr);
t2 =toc;

cout2 = interleaved2complex(coutPtr.value);
sizeaccum = 0;
for p=1:Psize
    nextsizeaccum = sizeaccum + M2(P(p))*N(P(p))*W;
    figure(4+p); plotdgtreal(reshape(postpad(cout2(1+sizeaccum:nextsizeaccum), M2(P(p))*N(P(p)) ),M2(P(p)),N(P(p))),a(P(p)),M(P(p)),'clim',[-90,10]);
    sizeaccum = nextsizeaccum;
    %ylim([0,0.01]);
end

% funname = makelibraryname('dgtrealmp_revert',flags.complexity,0);
% calllib('libltfat',funname,plan,coutPtr);

%%%%%%%%%%%%%%
errdb = libpointer('doublePtr',[1]);
funname = makelibraryname('dgtrealmp_get_errdb',flags.complexity,0);
calllib('libltfat',funname,plan,errdb);
err2 = errdb.value;
%%%%%%%%%%%%%%
%%%%%%%%%%%%%%
numitersPtr = libpointer('uint64Ptr',[1]);
funname = makelibraryname('dgtrealmp_get_numiters',flags.complexity,0);
calllib('libltfat',funname,plan,numitersPtr);
%%%%%%%%%%%%%%
%%%%%%%%%%%%%%
numatomsPtr = libpointer('uint64Ptr',[1]);
funname = makelibraryname('dgtrealmp_get_numatoms',flags.complexity,0);
calllib('libltfat',funname,plan,numatomsPtr);
%%%%%%%%%%%%%%

if doplot
cres2 = complex2interleaved(...
cast(randn(sizeaccum,1)+...
         1i*randn(sizeaccum,1),flags.complexity));
cresPtr = libpointer(dataPtr,cres2);

funname = makelibraryname('dgtrealmp_getresidualcoef_compact',flags.complexity,0);
calllib('libltfat',funname,plan,cresPtr);

cres = interleaved2complex(cresPtr.value);
sizeaccum = 0;
for p=1:Psize
    nextsizeaccum = sizeaccum + M2(P(p))*N(P(p))*W;
    figure(10+p); plotdgtreal(reshape(postpad(cres(1+sizeaccum:nextsizeaccum), M2(P(p))*N(P(p)) ),M2(P(p)),N(P(p))),1,100,'clim',[-90,10]);
    sizeaccum = nextsizeaccum;
end
end

%cres2 = reshape( postpad( interleaved2complex(cresPtr.value),M2(P(1))*N(P(1))  ),M2(P(1)),N(P(1)));
%figure(3); plotdgtreal(cres2,1,100,'clim',[-90,10]);
%figure(4); plotdgtreal(cres1-cres2,1,100,'dynrange',90);

fprintf('Init %.3f, reset %.3f, execute %.3f, both %.3f seconds, status %s.\n',tinit,t1,t2,t1+t2,dgtrealmpstring(statusExecute));

clear coutPtr cout
atoms = numel(find(abs(cout2(:))));

coutCell = cell(Psize,1);
sizeaccum = 0;
fout(:) = 0;
for p=1:Psize
    coutCell{p} = cout2(sizeaccum +1: sizeaccum + M2(P(p))*N(P(p))*W);
    coutCell{p} = reshape(coutCell{p},M2(P(p)),N(P(p)),W);
    sizeaccum = sizeaccum + M2(P(p))*N(P(p))*W;
    fout(:,p) = idgtreal(coutCell{p},gCell{p},a(P(p)),M(P(p)),mphaseconv);
end

indatoms = cellfun(@(a) numel(find(abs(a(:)))), coutCell);
atomstr = sprintf('%d,',indatoms);

clear cout2
figure(1);
plot((0:L-1)/fs,[fPtr.value, fout]);

errdb = 20*log10(norm(fPtr.value -sum(fout,2))/norm(fPtr.value));

errdb22 = libpointer(dataPtr,[1]);
foutPtr = libpointer(dataPtr,sum(fout,2));
funname = makelibraryname('snr',flags.complexity,0);
statusExecute = calllib('libltfat',funname,fPtr,foutPtr,L,errdb22);
%soundsc(fout,44100);pause(4);

errdb22.value;

fprintf('%i(%i) [%s] atoms (from %i), %i iters, sparsity %.3f, L=%i, %i atoms/s\nErr: True: %.8f dB, En: %.8f dB,\n',...
     atoms, numatomsPtr.value, atomstr(1:end-1),sum(M2(P(p))*L./a(P(p))),numitersPtr.value,atoms/L,L,atoms/t2, errdb22.value,err2);

shg;

funname = makelibraryname('dgtrealmp_done',flags.complexity,0);
statusDone = calllib('libltfat',funname,plan);

  atomsArr(end + 1) = atoms;
  [test_failed,fail]=ltfatdiditfail(abs(errdb - errtoldb)>0.1 ,test_failed);
    %[test_failed,fail]=ltfatdiditfail(res+statusInit,test_failed);
    %fprintf(['DGTREAL FREQINV WP auto %s L:%3i, W:%3i, a:%3i, M:%3i %s %s %s\n'],dirstr,L,W,a,M,flags.complexity,ltfatstatusstring(statusExecute),fail);
 %drawnow

end
end

function sstring=dgtrealmpstring(status)

[~,~,enuminfo]=libltfatprotofile;

map = structfun(@(a) a==status ,enuminfo.ltfat_dgtmp_status);
names = fieldnames(enuminfo.ltfat_dgtmp_status);
sstring = names{map};

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [bandwidth,density,X,Y]=kde2d(data,n,MIN_XY,MAX_XY)
% fast and accurate state-of-the-art
% bivariate kernel density estimator
% with diagonal bandwidth matrix.
% The kernel is assumed to be Gaussian.
% The two bandwidth parameters are
% chosen optimally without ever
% using/assuming a parametric model for the data or any "rules of thumb".
% Unlike many other procedures, this one
% is immune to accuracy failures in the estimation of
% multimodal densities with widely separated modes (see examples).
% INPUTS: data - an N by 2 array with continuous data
%            n - size of the n by n grid over which the density is computed
%                n has to be a power of 2, otherwise n=2^ceil(log2(n));
%                the default value is 2^8;
% MIN_XY,MAX_XY- limits of the bounding box over which the density is computed;
%                the format is:
%                MIN_XY=[lower_Xlim,lower_Ylim]
%                MAX_XY=[upper_Xlim,upper_Ylim].
%                The dafault limits are computed as:
%                MAX=max(data,[],1); MIN=min(data,[],1); Range=MAX-MIN;
%                MAX_XY=MAX+Range/4; MIN_XY=MIN-Range/4;
% OUTPUT: bandwidth - a row vector with the two optimal
%                     bandwidths for a bivaroate Gaussian kernel;
%                     the format is:
%                     bandwidth=[bandwidth_X, bandwidth_Y];
%          density  - an n by n matrix containing the density values over the n by n grid;
%                     density is not computed unless the function is asked for such an output;
%              X,Y  - the meshgrid over which the variable "density" has been computed;
%                     the intended usage is as follows:
%                     surf(X,Y,density)
% Example (simple Gaussian mixture)
% clear all
%   % generate a Gaussian mixture with distant modes
%   data=[randn(500,2);
%       randn(500,1)+3.5, randn(500,1);];
%   % call the routine
%     [bandwidth,density,X,Y]=kde2d(data);
%   % plot the data and the density estimate
%     contour3(X,Y,density,50), hold on
%     plot(data(:,1),data(:,2),'r.','MarkerSize',5)
%
% Example (Gaussian mixture with distant modes):
%
% clear all
%  % generate a Gaussian mixture with distant modes
%  data=[randn(100,1), randn(100,1)/4;
%      randn(100,1)+18, randn(100,1);
%      randn(100,1)+15, randn(100,1)/2-18;];
%  % call the routine
%    [bandwidth,density,X,Y]=kde2d(data);
%  % plot the data and the density estimate
%  surf(X,Y,density,'LineStyle','none'), view([0,60])
%  colormap hot, hold on, alpha(.8)
%  set(gca, 'color', 'blue');
%  plot(data(:,1),data(:,2),'w.','MarkerSize',5)
%
% Example (Sinusoidal density):
%
% clear all
%   X=rand(1000,1); Y=sin(X*10*pi)+randn(size(X))/3; data=[X,Y];
%  % apply routine
%  [bandwidth,density,X,Y]=kde2d(data);
%  % plot the data and the density estimate
%  surf(X,Y,density,'LineStyle','none'), view([0,70])
%  colormap hot, hold on, alpha(.8)
%  set(gca, 'color', 'blue');
%  plot(data(:,1),data(:,2),'w.','MarkerSize',5)
%
%  Reference:
% Kernel density estimation via diffusion
% Z. I. Botev, J. F. Grotowski, and D. P. Kroese (2010)
% Annals of Statistics, Volume 38, Number 5, pages 2916-2957.
global N A2 I
if nargin<2
    n=2^8;
end
n=2^ceil(log2(n)); % round up n to the next power of 2;
N=size(data,1);
if nargin<3
    MAX=max(data,[],1); MIN=min(data,[],1); Range=MAX-MIN;
    MAX_XY=MAX+Range/2; MIN_XY=MIN-Range/2;
end
scaling=MAX_XY-MIN_XY;
if N<=size(data,2)
    error('data has to be an N by 2 array where each row represents a two dimensional observation')
end
transformed_data=(data-repmat(MIN_XY,N,1))./repmat(scaling,N,1);
%bin the data uniformly using regular grid;
initial_data=ndhist(transformed_data,n);
% discrete cosine transform of initial data
a= dct2d(initial_data);
% now compute the optimal bandwidth^2
  I=(0:n-1).^2; A2=a.^2;
 t_star=root(@(t)(t-evolve(t)),N);
p_02=func([0,2],t_star);p_20=func([2,0],t_star); p_11=func([1,1],t_star);
t_y=(p_02^(3/4)/(4*pi*N*p_20^(3/4)*(p_11+sqrt(p_20*p_02))))^(1/3);
t_x=(p_20^(3/4)/(4*pi*N*p_02^(3/4)*(p_11+sqrt(p_20*p_02))))^(1/3);
% smooth the discrete cosine transform of initial data using t_star
a_t=exp(-(0:n-1)'.^2*pi^2*t_x/2)*exp(-(0:n-1).^2*pi^2*t_y/2).*a;
% now apply the inverse discrete cosine transform
if nargout>1
    density=idct2d(a_t)*(numel(a_t)/prod(scaling));
	density(density<0)=eps; % remove any negative density values
    [X,Y]=meshgrid(MIN_XY(1):scaling(1)/(n-1):MAX_XY(1),MIN_XY(2):scaling(2)/(n-1):MAX_XY(2));
end
bandwidth=sqrt([t_x,t_y]).*scaling;
end
%#######################################
function  [out,time]=evolve(t)
global N
Sum_func = func([0,2],t) + func([2,0],t) + 2*func([1,1],t);
time=(2*pi*N*Sum_func)^(-1/3);
out=(t-time)/time;
end
%#######################################
function out=func(s,t)
global N
if sum(s)<=4
    Sum_func=func([s(1)+1,s(2)],t)+func([s(1),s(2)+1],t); const=(1+1/2^(sum(s)+1))/3;
    time=(-2*const*K(s(1))*K(s(2))/N/Sum_func)^(1/(2+sum(s)));
    out=psi(s,time);
else
    out=psi(s,t);
end

end
%#######################################
function out=psi(s,Time)
global I A2
% s is a vector
w=exp(-I*pi^2*Time).*[1,.5*ones(1,length(I)-1)];
wx=w.*(I.^s(1));
wy=w.*(I.^s(2));
out=(-1)^sum(s)*(wy*A2*wx')*pi^(2*sum(s));
end
%#######################################
function out=K(s)
out=(-1)^s*prod((1:2:2*s-1))/sqrt(2*pi);
end
%#######################################
function data=dct2d(data)
% computes the 2 dimensional discrete cosine transform of data
% data is an nd cube
[nrows,ncols]= size(data);
if nrows~=ncols
    error('data is not a square array!')
end
% Compute weights to multiply DFT coefficients
w = [1;2*(exp(-i*(1:nrows-1)*pi/(2*nrows))).'];
weight=w(:,ones(1,ncols));
data=dct1d(dct1d(data)')';
    function transform1d=dct1d(x)

        % Re-order the elements of the columns of x
        x = [ x(1:2:end,:); x(end:-2:2,:) ];

        % Multiply FFT by weights:
        transform1d = real(weight.* fft(x));
    end
end
%#######################################
function data = idct2d(data)
% computes the 2 dimensional inverse discrete cosine transform
[nrows,ncols]=size(data);
% Compute wieghts
w = exp(i*(0:nrows-1)*pi/(2*nrows)).';
weights=w(:,ones(1,ncols));
data=idct1d(idct1d(data)');
    function out=idct1d(x)
        y = real(ifft(weights.*x));
        out = zeros(nrows,ncols);
        out(1:2:nrows,:) = y(1:nrows/2,:);
        out(2:2:nrows,:) = y(nrows:-1:nrows/2+1,:);
    end
end
%#######################################
function binned_data=ndhist(data,M)
% this function computes the histogram
% of an n-dimensional data set;
% 'data' is nrows by n columns
% M is the number of bins used in each dimension
% so that 'binned_data' is a hypercube with
% size length equal to M;
[nrows,ncols]=size(data);
bins=zeros(nrows,ncols);
for i=1:ncols
    [dum,bins(:,i)] = histc(data(:,i),[0:1/M:1],1);
    bins(:,i) = min(bins(:,i),M);
end
% Combine the  vectors of 1D bin counts into a grid of nD bin
% counts.
binned_data = accumarray(bins(all(bins>0,2),:),1/nrows,M(ones(1,ncols)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function t=root(f,N)
% try to find smallest root whenever there is more than one
N=50*(N<=50)+1050*(N>=1050)+N*((N<1050)&(N>50));
tol=10^-12+0.01*(N-50)/1000;
flag=0;
while flag==0
    try
        t=fzero(f,[0,tol]);
        flag=1;
    catch
        tol=min(tol*2,.1); % double search interval
    end
    if tol==.1 % if all else fails
        t=fminbnd(@(x)abs(f(x)),0,.1); flag=1;
    end
end
end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
%kmscalar_sub 'Scalar Subtract: Output = Input Matrix - Constant'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros mscalar_sub.pane file
%
% Parameters:
% InputFile: i1 'Input Matrix', required: 'Input matrix'
% OutputFile: o 'Output ', required: 'Output matrix'
%
% Example: o = kmscalar_sub(i1, {'i1','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% mscalar_sub - Output = Input Matrix - Constant
%
%  DESCRIPTION
% .I mscalar_sub
% subtracts a specified constant from each component in a matrix.
% "mscalar_sub\fR is really a pane object that uses the \fIkarith2\fR operator
% in the DATAMANIP toolbox to actually perform the operation, so the man page for
% "karith2\fR should be consulted for additional details regarding the treatment
% of maps, masks, etc.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% DATAMANIP::karith2(1)  DATAMANIP::ksub(1)
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kmscalar_sub(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kmscalar_sub(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'o', '__output'};
maxval={0,0};
minval={0,0};
istoggle=[0,0];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'karith2"  -sub'],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
clear all

nInst = 500;
nVars = 200;
X = randn(nInst,nVars);
w = randn(nVars,1);
y = sign(X*w + randn(nInst,1));

w_init = zeros(nVars,1);
funObj = @(w)LogisticLoss(w,X,y);

fprintf('\nRunning Steepest Descent\n');
options.Method = 'sd';
minFunc(@LogisticLoss,w_init,options,X,y);
pause;

fprintf('\nRunning Cyclic Steepest Descent\n');
options.Method = 'csd';
minFunc(@LogisticLoss,w_init,options,X,y);
pause;

fprintf('\nRunning Conjugate Gradient\n');
options.Method = 'cg';
minFunc(@LogisticLoss,w_init,options,X,y);
pause;

fprintf('\nRunning Scaled Conjugate Gradient\n');
options.Method = 'scg';
minFunc(@LogisticLoss,w_init,options,X,y);
pause;

fprintf('\nRunning Preconditioned Conjugate Gradient (Diagonal preconditioner)\n');
options.Method = 'pcg';
options.precFunc = @LogisticDiagPrecond;
minFunc(@LogisticLoss,w_init,options,X,y);
pause;

fprintf('\nRunning Preconditioned Conjugate Gradient (L-BFGS preconditioner)\n');
options.Method = 'pcg';
options.precFunc = [];
minFunc(@LogisticLoss,w_init,options,X,y);
pause;

fprintf('\nRunning Hessian-Free Newton w/ numerical Hessian-Vector products\n');
options.Method = 'newton0';
minFunc(@LogisticLoss,w_init,options,X,y);
pause;

fprintf('\nRunning Preconditioned Hessian-Free Newton w/ numerical Hessian-Vector products (Diagonal preconditioner)\n');
options.Method = 'pnewton0';
options.precFunc = @LogisticDiagPrecond;
minFunc(@LogisticLoss,w_init,options,X,y);
pause;

fprintf('\nRunning Preconditioned Hessian-Free Newton w/ numerical Hessian-Vector products (L-BFGS preconditioner)\n');
options.Method = 'pnewton0';
options.precFunc = [];
minFunc(@LogisticLoss,w_init,options,X,y);
pause;

fprintf('\nRunning Hessian-Free Newton w/ analytic Hessian-Vector products\n');
options.Method = 'newton0';
options.HvFunc = @LogisticHv;
minFunc(@LogisticLoss,w_init,options,X,y);
pause;

fprintf('\nRunning Preconditioned Hessian-Free Newton w/ analytic Hessian-Vector products (Diagonal preconditioner)\n');
options.Method = 'pnewton0';
options.HvFunc = @LogisticHv;
options.precFunc = @LogisticDiagPrecond;
minFunc(@LogisticLoss,w_init,options,X,y);
pause;

fprintf('\nRunning Preconditioned Hessian-Free Newton w/ analytic Hessian-Vector products (L-BFGS preconditioner)\n');
options.Method = 'pnewton0';
options.precFunc = [];
options.HvFunc = @LogisticHv;
minFunc(@LogisticLoss,w_init,options,X,y);
pause;

author: zouchuhang


text:
function x=eqConstLSSpher(A,b,alpha,epsRed)
%%EQCONSTLSSPHER Find x to minimize norm(A*x-b,2) under the constraint
%                that norm(x,2)=alpha. This is essentialy contraining x to
%                the surface of a sphere of radius alpha. This only solves
%                real systems.
%
%INPUTS: A A real mXn matrix with m>=n.
%        b A real mX1 vector.
%    alpha The equality constraint value. If this parameter is omitted or
%          an empty matrix is passed, the default of 1 is used.
%   epsRed Let xU be the unconstrained solution to the problem. A possible
%          constrained solutions is considered valid if
%          abs(norm(x)^2-alpha^2)<=epsRed*abs(norm(xU)^2-alpha^2)
%          If no such solutions are found, then xU is returned. The default
%          for this parameter if omitted or an empty matrix is passed is
%          1e-9. This parameter should be between 0 and 1.
%
%OUTPUTS: x The optimal value of x subject to the spherical constraint.
%    lambda The Lagrangian multiplier used in the optimization. lambda>=0.
%           If lambda=0, then the constraint on x did not have to be
%           enforced.
%
%This implements a modified version of the algorithm of Chapter 6.2.1 of
%[1]. The algorithm of Chapter 6.2.1 of [1] solves the optimization
%constrained such that norm(x,2)<=alpha. We wish to solve the equality
%constrained problem. To do so, we always enforce the constraint. An
%equation has to be solved for scalar zeros over lambda. To do so, we
%multiply both sids by the denominators resulting in a polynomial system.
%The system is then solved. However, some solutions might not satisfy the
%constraint (due to cancellation in denominators). Thus, candidate
%solutions that do not improve the error in the magnitude by a sufficient
%amount compared to the unconstrained solution are discarded. If no
%solutions are left, the the unconstrained solution is used.
%
%EXAMPLE:
%This is a simple example where the x returned by the
%inequality-constrained algorithm is too small.
% A=magic(8)+20*eye(8);
% b=(1:8).';
% x=inv(A)*b;
% norm(x)%Norm <1.
% xConst=eqConstLSSpher(A,b);
% norm(xConst)%Norm=1
% %It is not just normalizing the vector; elements are not scaled by a
% %constant.
% x./xConst
%
%REFERENCES:
%[1] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed.
%    Baltimore: Johns Hopkins University Press, 2013.
%
%December 2020 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(alpha))
    alpha=1;
end

if(nargin<4||isempty(epsRed))
    epsRed=1e-9;
end

r=rank(A);

[U,Sigma,V]=svd(A,0);
sigma=diag(Sigma);

%Sums are only up to r, so get rid of the extra elements.
U=U(:,1:r);
V=V(:,1:r);
sigma=sigma(1:r);

bTilde=U'*b;
xUnconst=V*(bTilde./sigma);

alpha2Unconst=abs(norm(xUnconst)^2-alpha^2);

numVal=(sigma.*bTilde).^2;
denomVal=sigma.^2;

numDim=length(sigma);
polynoms=[ones(numDim,1),2*denomVal,denomVal.^2];

%Construct the polynomial to solve.
polyNom=0;
for k1=1:numDim
    curPoly=numVal(k1);
    for k2=1:numDim
        if(k2==k1)
            continue;
        end
        curPoly=conv(curPoly,polynoms(k2,:));
    end
    polyNom=polySum(curPoly,polyNom);
end

%The final term
curPoly=-alpha^2;
for k2=1:numDim
    curPoly=conv(curPoly,polynoms(k2,:));
end
polyNom=polySum(curPoly,polyNom);
lambdaVals=roots(polyNom).';
%Get rid of imaginary solutions.
lambdaVals=lambdaVals(imag(lambdaVals)==0);
numSol=length(lambdaVals);

%Due to certain values corresponding to zero denominators, there can be
%more solutions in lambda than are valid. We must eliminate all solutions
%that do not produce x values with the correct magnitude. First, get all
%possibly valid solutions regardless of the magnitude.
x=zeros(numDim,numSol);
numKept=0;
for k=1:numSol
    xCur=V*((sigma.*bTilde)./(sigma.^2+lambdaVals(k)));
    normErr=abs(norm(xCur)^2-alpha^2);
    if(all(isfinite(xCur))&&normErr<=epsRed*alpha2Unconst)
        numKept=numKept+1;
        x(:,numKept)=xCur;
    end
end

if(numKept==0)
    %If nothing was kept, then just use the unconstrained solution.
    x=xUnconst;
    return
end
x=x(:,1:numKept);

%Of all of the solutions kept, take the one that minimizes the original
%optimization problem.
if(numKept>1)
    minCost=norm(A*x(:,1)-b,2);
    minIdx=1;
    for k=2:numKept
       curCost=norm(A*x(:,k)-b,2);

       if(curCost<minCost)
           minIdx=k;
           minCost=curCost;
       end
    end
    x=x(:,minIdx);
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 plotDisturbanceComparison(files, laps, labels, max_data_points)
% files: .mat files in a dict. like [controller1 datapoint 1, controller1
% datapoint 2, ...., controller2 datapoint 1,...., controllerN datapoint
% max_data_points]
% max_data_points: number of data points per controller

LineWidth = 1;

k = length(files)/max_data_points;
lap_times_flat = zeros(length(files), 1);
acc_viol_flat = zeros(length(files), 1);
dy_viol_flat = zeros(length(files), 1);
sum_viol_flat = zeros(length(files), 1);
dax_flat = zeros(length(files), 1);
lap_times = zeros(k, max_data_points);
acc_viol = zeros(k, max_data_points);
dy_viol = zeros(k, max_data_points);
sum_viol = zeros(k, max_data_points);
dax = zeros(k,max_data_points);

for i = 1:1:length(files)
    data = load(files{i});
    idx_start = find((data.debug.debug_mvdc_path_matching_debug_ActualTrajPoint_LapCnt.Data == laps(1)), 1, 'first');
    idx_end = find((data.debug.debug_mvdc_path_matching_debug_ActualTrajPoint_LapCnt.Data == laps(2)), 1, 'last');
    if(isempty(idx_end))
        disp(['No lap found for dataset ' files{i}]);
        lap_times_flat(i) = NaN;
        acc_viol_flat(i) = NaN;
        dy_viol_flat(i) = NaN;
        dax_flat(i) = NaN;
    else
        lap_times_flat(i) = data.debug.debug_mvdc_trajectory_driver_perf_LapTime_s.Data(idx_end);
        acc_viol_flat(i) = data.debug.debug_mvdc_trajectory_driver_perf_TireConstViolations.Data(idx_end);
        dy_viol_flat(i) = data.debug.debug_mvdc_trajectory_driver_perf_LatConstViolations.Data(idx_end);
        dax_flat(i) = data.debug_slow.debug_slow_M_tilde.Data(500,1);
    end
end

% resort to generate appropriate plots
for i = 1:1:k
    for j = 1:1:max_data_points
        lap_times(i,j) = lap_times_flat((i-1)*max_data_points+j);
        acc_viol(i,j) = acc_viol_flat((i-1)*max_data_points+j);
        dy_viol(i,j) = dy_viol_flat((i-1)*max_data_points+j);
        dax(i,j) = dax_flat((i-1)*max_data_points+j);
    end
end

figure;
subplot(3, 1, 1); box on, hold on, grid on;
for i = 1:1:k
    plot(dax(i,:), dy_viol(i,:),'LineWidth', LineWidth);
end
xlabel('dax in mps2'); ylabel('violated dy constraints');
legend(labels{1:end});

subplot(3, 1, 2); box on, hold on, grid on;
for i = 1:1:k
    plot(dax(i,:), acc_viol(i,:),'LineWidth', LineWidth);
end
xlabel('dax in mps2'); ylabel('violated a constraints');

subplot(3, 1, 3); box on, hold on, grid on;
for i = 1:1:k
    plot(dax(i,:), lap_times(i,:),'LineWidth', LineWidth);
end
xlabel('dax in mps2'); ylabel('Lap time in s');

author: TUMFTM


text:
% [PYR, INDICES, STEERMTX, HARMONICS] = buildSpyr(IM, HEIGHT, FILTFILE, EDGES)
%
% Construct a steerable pyramid on matrix IM.  Convolutions are
% done with spatial filters.
%
% HEIGHT (optional) specifies the number of pyramid levels to build. Default
% is maxPyrHt(size(IM),size(FILT)).
% You can also specify 'auto' to use this value.
%
% FILTFILE (optional) should be a string referring to an m-file that
% returns the rfilters.  (examples: 'sp0Filters', 'sp1Filters',
% 'sp3Filters','sp5Filters'.  default = 'sp1Filters'). EDGES specifies
% edge-handling, and defaults to 'reflect1' (see corrDn).
%
% PYR is a vector containing the N pyramid subbands, ordered from fine
% to coarse.  INDICES is an Nx2 matrix containing the sizes of
% each subband.  This is compatible with the MatLab Wavelet toolbox.
% See the function STEER for a description of STEERMTX and HARMONICS.

% Eero Simoncelli, 6/96.
% See http://www.cis.upenn.edu/~eero/steerpyr.html for more
% information about the Steerable Pyramid image decomposition.

function [pyr,pind,steermtx,harmonics] = buildSpyr(im, ht, filtfile, edges)

%-----------------------------------------------------------------
%% DEFAULTS:

if (exist('filtfile') ~= 1)
  filtfile = 'sp1Filters';
end

if (exist('edges') ~= 1)
  edges= 'reflect1';
end

if (isstr(filtfile) & (exist(filtfile) == 2))
   [lo0filt,hi0filt,lofilt,bfilts,steermtx,harmonics] = eval(filtfile);
else
  fprintf(1,'\nUse buildSFpyr for pyramids with arbitrary numbers of orientation bands.\n');
  error('FILTFILE argument must be the name of an M-file containing SPYR filters.');
end

max_ht = maxPyrHt(size(im), size(lofilt,1));
if ( (exist('ht') ~= 1) | (ht == 'auto') )
  ht = max_ht;
else
  if (ht > max_ht)
    error(sprintf('Cannot build pyramid higher than %d levels.',max_ht));
  end
end

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

hi0 = corrDn(im, hi0filt, edges);
lo0 = corrDn(im, lo0filt, edges);

[pyr,pind] = buildSpyrLevs(lo0, ht, lofilt, bfilts, edges);

pyr = [hi0(:) ; pyr];
pind = [size(hi0); pind];

repo: SelfExSR


text:
function [M,Id] = compute_texture_patchwork(H,n, options)

% compute_texture_patchwork - mix several textures
%
%   M = compute_texture_patchwork(H,n,options);
%
%   H can be a cell array of textures or a (n,n,3) matrix.
%   Handles up to 5 textures.
%
%   To generate a random partition, set options.patchwork_mode='random'.
%
%   Copyright (c) 2006 Gabriel Peyre

if not(iscell(H))
    H1 = H; H = {};
    for i=1:size(H1,3)
        H{i} = H1(:,:,i);
    end
end

options.null = 0;
mode = getoptions(options, 'patchwork_mode', 'deterministic');

if nargin<2
    n = size(H{1},1);
end

p = length(H); % number of textures

if strcmp(mode, 'random')

    sigma = getoptions(options, 'patchwork_sigma', 60);
    J = floor(log2(p));
    A = ones(n);
    for j=0:J-1
        B = A;
        for k=1:2^j
            I = find(B==k);
            U = perform_blurring(randn(n),sigma,options);
            s = median(U(I));
            I1 = find( (B==k) & (U>s) );
            I2 = find( (B==k) & (U<=s) );
            A(I1) = 2*k-1;
            A(I2) = 2*k;
        end
    end

    M = zeros(n);
    Id = zeros(n);
    for i=1:p
        B = H{i};
        if size(B,1)<n
            B = perform_image_extension(B,n);
        end
        B = B(1:n,1:n);
        M(A==i) = B(A==i);
        Id(A==i) = i;
    end
    return;

end

m = min(size(H{1},1),n);

M = zeros(n);
M(1:m,1:m) = H{1}(1:m,1:m);
Id = ones(n);

if p==2
    M(end/2+1:end,:) = H{2}(end-n/2+1:end,:);
    Id(end/2+1:end,:) = 2;
    return;
end

if p>1
    M(1:end/2,end/2+1:end) = H{2}(1:1:n/2,1:n/2);
    Id(1:end/2,end/2+1:end) = 2;
end
if p>2
    M(end/2+1:end,1:end/2) = H{3}(1:1:n/2,1:n/2);
    Id(end/2+1:end,1:end/2) = 3;
end
if p>3
    M(end/2+1:end,end/2+1:end) = H{4}(1:1:n/2,1:n/2);
    Id(end/2+1:end,end/2+1:end) = 4;
end
if p>4
    r = 1;
    A = H{5}(1:1:n/2,1:n/2);
    x = linspace(-1,1,n/2)';
    [Y,X] = meshgrid(x,x);
    J = find( X.^2 + Y.^2 <= r^2 );
    x  = [ones(n/4,1)*Inf; x; ones(n/4,1)*Inf];
    [Y,X] = meshgrid(x,x);
    I = find( X.^2 + Y.^2 <= r^2 );
    M(I) = A(J);
    Id(I) = 5;
end

author: gpeyre


text:
function x = RofCurve(X,Y)

% X=[10 0 -8];
% Y=[0 10 0];

% Defining the three nonlinear equations to be solved (equation of circle
% passing through the 3 points)
defn1=['[ ( x(1)-' num2str(X(1)) ' )^2+( x(2)-' num2str(Y(1)) ')^2-x(3)^2 ;' ];
defn2=[  '( x(1)-' num2str(X(2)) ' )^2+( x(2)-' num2str(Y(2)) ')^2-x(3)^2 ;' ];
defn3=[  '( x(1)-' num2str(X(3)) ' )^2+( x(2)-' num2str(Y(3)) ')^2-x(3)^2 ]' ];
inlinedef=[defn1 defn2 defn3];

%Defining the inline function that will compute  the error, which has to be
%minimized
myfun=inline(inlinedef,'x');

options = optimset('Display','off','TolX',1e-10);
x0=[0;0;0];
[x,fval] = fsolve(myfun,x0,options);
x=[x(1) x(2) abs(x(3))];

author: Sable

repo: mcbench-benchmarks