text:
function calc_bootfitloglike_a2(a,time,timef,bootloops,maepi)
    % function calc_bootfitloglike_a2(a,time,timef,bootloops,maepi);
    % --------------------------------------------------
    % Plots Ncum observed vs. Ncum modeled for specified time windows
    %
    % Input variables:
    % a         : earthquake catalog
    % time      : learning period fo fit Omori parameters
    % timef     : forecast period
    % bootloops : Number of bootstraps
    % maepi     : mainshock
    %
    % J.Woessner, S. Wiemer
    % last update: 05.08.03

report_this_filefun(mfilename('fullpath'));
    % Surpress warnings from fmincon
    warning off;

report_this_filefun(mfilename('fullpath'));
    %[m_main, main] = max(a.Magnitude);
    date_matlab = datenum(a.Date.Year,a.Date.Month,a.Date.Day,a.Date.Hour,a.Date.Minute,zeros(size(a,1),1));
    date_main = datenum(floor(maepi(3)),maepi(4),maepi(5),maepi(8),maepi(9),0);
    time_aftershock = date_matlab-date_main;
    % Select biggest aftershock earliest in time, but more than 1 day after mainshock
    fDay = 1;
    ft_c=fDay/365; % Time not considered to find biggest aftershock
    vSel = (a.Date > maepi(:,3)+ft_c & a.Date<= maepi(:,3)+time/365);
    mCat = a.subset(vSel);
    vSel = mCat(:,6) == max(mCat(:,6));
    vBigAf = mCat(vSel,:);
    if length(mCat(:,1)) > 1
        vSel = vBigAf(:,3) == min(vBigAf(:,3));
        vBigAf = vBigAf(vSel,:);
    end

    date_biga = datenum(floor(vBigAf(3)),vBigAf(4),vBigAf(5),vBigAf(8),vBigAf(9),0);
    fT1 = date_biga - date_main; % Time of big aftershock

    % Aftershock times
    l = time_aftershock(:) > 0;
    tas = time_aftershock(l);
    eqcatalogue = a.subset(l);

    % time_as: Learning period
    l = tas <= time;
    time_as=tas(l);

    % Times up to the forecast time
    lf = tas <= time+timef ;
    time_asf= [tas(lf) ];
    time_asf=sort(time_asf);

    % Calculate p,c,k for dataset
    prompt  = {'Enter model number (1:pck, 2:pckk, 3:ppckk, 4:ppcckk:'};
    title   = 'Model selection for fitting aftershock sequence';
    lines= 1;
    def     = {'1'};
    answer  = inputdlg(prompt,title,lines,def);
    nMod = str2double(answer{1});

    % Calculate uncertainty and mean values of p,c,and k
    [mMedModF, mStdL, loopout] = brutebootloglike_a2(time_as, time_asf, bootloops,fT1,nMod);
    pmed1 = mMedModF(1,1);
    pmed2 = mMedModF(1,3);
    cmed1 = mMedModF(1,5);
    cmed2 = mMedModF(1,7);
    kmed1 = mMedModF(1,9);
    kmed2 = mMedModF(1,11);

    % Compute model according to model choice
    if nMod == 1
        [pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL] = bruteforceloglike_a2(time_as,fT1,nMod);
    elseif nMod == 2
        [pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL] = bruteforceloglike_a2(time_as,fT1,nMod);
    elseif nMod == 3
        [pval1, pval2, cval1, cval2, kval1, kval2 , fAIC, fL] = bruteforceloglike_a2(time_as,fT1,nMod);
    else
        [pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL] = bruteforceloglike_a2(time_as,fT1,nMod);
    end
    % Round values
    pval1 = round(100*pval1)/100;
    pval2 = round(100*pval2)/100;
    cval1 = round(100*cval1)/100;
    cval2 = round(100*cval2)/100;
    kval1 = round(10*kval1)/10;
    kval2 = round(10*kval2)/10;

    if (isnan(pval1) == 0 & isnan(pval2) == 0)

        figure_w_normalized_uicontrolunits('Numbertitle','off','Name','Forecast aftershock occurence')
        loopout = [loopout , loopout(:,1)*0];

        % Time until end of forecast
        for j = 1:length(loopout(:,1))
            cumnr = (1:length(time_asf))';
            cumnr_model = [];
            pval1 = loopout(j,1);
            pval2 = loopout(j,2);
            cval1 = loopout(j,3);
            cval2 = loopout(j,4);
            kval1 = loopout(j,5);
            kval2 = loopout(j,6);
            if nMod == 1
                for i=1:length(time_asf)
                    if pval1 ~= 1
                        cm = kval1/(pval1-1)*(cval1^(1-pval1)-(time_asf(i)+cval1)^(1-pval1));
                    else
                        cm = kval1*log(time_asf(i)/cval1+1);
                    end
                    cumnr_model = [cumnr_model; cm];
                end % END of FOR on length(time_asf)
                loopout(j,9) = max(cumnr_model);
            else
                for i=1:length(time_asf)
                    if time_asf(i) <= fT1
                        if pval1 ~= 1
                            cm = kval1/(pval1-1)*(cval1^(1-pval1)-(time_asf(i)+cval1)^(1-pval1));
                        else
                            cm = kval1*log(time_asf(i)/cval1+1);
                        end
                        cumnr_model = [cumnr_model; cm];
                    else
                        if (pval1 ~= 1 & pval2 ~= 1)
                            cm = kval1/(pval1-1)*(cval1^(1-pval1)-(time_asf(i)+cval1)^(1-pval1))+ kval2/(pval2-1)*(cval2^(1-pval2)-(time_asf(i)-fT1+cval2)^(1-pval2));
                        else
                            cm = kval1*log(time_asf(i)/cval1+1) + kval2*log((time_asf(i)-fT1)/cval2+1);
                        end
                        cumnr_model = [cumnr_model; cm];
                    end %END of IF on fT1
                end % End of FOR length(time_asf)
                loopout(j,9) = max(cumnr_model);
            end % End of if on nMod
            pfloop = plot(time_asf,cumnr_model,'color',[0.8 0.8 0.8]);
            hold on
            %drawnow
        end
        % 2nd moment of bootstrap number of forecasted number of events
        fStdBst = calc_StdDev(loopout(:,9));
        %
        % Plot the forecast ...
        cumnrf = (1:length(time_asf))';
        cumnr_modelf = [];
        if nMod == 1
            for i=1:length(time_asf)
                if pval1 ~= 1
                    cm = kval1/(pval1-1)*(cval1^(1-pval1)-(time_asf(i)+cval1)^(1-pval1));
                else
                    cm = kval1*log(time_asf(i)/cval1+1);
                end
                cumnr_modelf = [cumnr_modelf; cm];
            end % END of FOR on length(time_asf)
        else
            for i=1:length(time_asf)
                if time_asf(i) <= fT1
                    if pval1 ~= 1
                        cm = kval1/(pval1-1)*(cval1^(1-pval1)-(time_asf(i)+cval1)^(1-pval1));
                    else
                        cm = kval1*log(time_asf(i)/cval1+1);
                    end
                    cumnr_modelf = [cumnr_modelf; cm];
                else
                    if (pval1 ~= 1 & pval2 ~= 1)
                        cm = kval1/(pval1-1)*(cval1^(1-pval1)-(time_asf(i)+cval1)^(1-pval1))+ kval2/(pval2-1)*(cval2^(1-pval2)-(time_asf(i)-fT1+cval2)^(1-pval2));
                    else
                        cm = kval1*log(time_asf(i)/cval1+1) + kval2*log((time_asf(i)-fT1)/cval2+1);
                    end
                    cumnr_modelf = [cumnr_modelf; cm];
                end %END of IF on fT1
            end % End of FOR length(time_asf)
        end % End of if on nMod
        time_asf=sort(time_asf);
        cumnr_modelf=sort(cumnr_modelf);

        pf1 =  plot(time_asf,cumnr_modelf,'g-.','Linewidth',2);
        hold on
        %pf2 =  plot(time_asf,cumnrf, 'b-','Linewidth',2);
        %
        % Plot the  fit to the observed data
        % Cumulative number of observed events
        cumnr = (1:length(time_as))';
        cumnr_model = [];
        if nMod == 1
            for i=1:length(time_as)
                if pval1 ~= 1
                    cm = kval1/(pval1-1)*(cval1^(1-pval1)-(time_as(i)+cval1)^(1-pval1));
                else
                    cm = kval1*log(time_as(i)/cval1+1);
                end
                cumnr_model = [cumnr_model; cm];
            end % END of FOR on length(time_as)
        else
            for i=1:length(time_as)
                if time_as(i) <= fT1
                    if pval1 ~= 1
                        cm = kval1/(pval1-1)*(cval1^(1-pval1)-(time_as(i)+cval1)^(1-pval1));
                    else
                        cm = kval1*log(time_as(i)/cval1+1);
                    end
                    cumnr_model = [cumnr_model; cm];
                else
                    if (pval1 ~= 1 & pval2 ~= 1)
                        cm = kval1/(pval1-1)*(cval1^(1-pval1)-(time_as(i)+cval1)^(1-pval1))+ kval2/(pval2-1)*(cval2^(1-pval2)-(time_as(i)-fT1+cval2)^(1-pval2));
                    else
                        cm = kval1*log(time_as(i)/cval1+1) + kval2*log((time_as(i)-fT1)/cval2+1);
                    end
                    cumnr_model = [cumnr_model; cm];
                end %END of IF on fT1
            end % End of FOR length(time_as)
        end % End of if on nMod
        time_as=sort(time_as);
        cumnr_model=sort(cumnr_model);
        p1 = plot(time_as,cumnr_model,'r','Linewidth',2,'Linestyle','--');
        hold on;
        p2 = plot(time_as,cumnr,'b','Linewidth',2,'Linestyle','--');

        % Plot the forecast from median value
        cumnr_modelmed = [];
        if nMod == 1
            for i=1:length(time_asf)
                if pval1 ~= 1
                    cm = kmed1/(pmed1-1)*(cmed1^(1-pmed1)-(time_asf(i)+cmed1)^(1-pmed1));
                else
                    cm = kmed1*log(time_asf(i)/cmed1+1);
                end
                cumnr_modelmed = [cumnr_modelmed; cm];
            end % END of FOR on length(time_asf)
        else
            for i=1:length(time_asf)
                if time_asf(i) <= fT1
                    if pmed1 ~= 1
                        cm = kmed1/(pmed1-1)*(cmed1^(1-pmed1)-(time_asf(i)+cmed1)^(1-pmed1));
                    else
                        cm = kmed1*log(time_asf(i)/cmed1+1);
                    end
                    cumnr_modelmed = [cumnr_modelmed; cm];
                else
                    if (pmed1 ~= 1 & pmed2 ~= 1)
                        cm = kmed1/(pmed1-1)*(cmed1^(1-pmed1)-(time_asf(i)+cmed1)^(1-pmed1))+ kmed2/(pmed2-1)*(cmed2^(1-pmed2)-(time_asf(i)-fT1+cmed2)^(1-pmed2));
                    else
                        cm = kmed1*log(time_asf(i)/cmed1+1) + kmed2*log((time_asf(i)-fT1)/cmed2+1);
                    end
                    cumnr_modelmed = [cumnr_modelmed; cm];
                end %END of IF on fT1
            end % End of FOR length(time_asf)
        end % End of if on nMod
        time_asf=sort(time_asf);
        cumnr_modelmed=sort(cumnr_modelmed);
        pmedmod =  plot(time_asf,cumnr_modelmed,'y-.','Linewidth',2);

        % Plot observed events in forecast period from endpoint of modeled events in learning period
        vSel = time_asf >= max(time_as);
        vCumnr_forecast = cumnrf(vSel,:);
        vTime_forecast = time_asf(vSel,:);
        % Difference of modelled and observed number of events at time_as
        fDiff_timeas = cumnr_modelmed(length(time_as))-cumnrf(length(time_as));
        vCumnr_forecast = vCumnr_forecast+fDiff_timeas;
        pf3 = plot(vTime_forecast, vCumnr_forecast,'m-.','Linewidth',2);

        xlabel('Time [days]')
        ylabel('Cumulative number of aftershocks')
        xlim([0 max(time_asf)]);

        % Plot standard deviation from bootstrap
        ps2=errorbar(max(time_asf),max(cumnr_modelmed),fStdBst,fStdBst);
        set(ps2,'Linewidth',2,'Color',[1 0 0])
        %
        legend([p2 p1 pf1 pf3 pmedmod min(ps2)],'data','model to data','forecast','observed', 'Median Bst-model', '\sigma (Bst)','location','Best');
        %
        %     % Title
        % Find amount of events in forecast period for modeled data
        nummod2 = max(cumnr_modelf)-cumnr_modelf(length(time_as));
        nummod = max(cumnr_modelmed)-cumnr_modelmed(length(time_as));
        %     % Find amount of  events in forecast period for observed data
        l = time_asf <=time+timef & time_asf > time;
        numreal = sum(l); % observed number of aftershocks
        %     fRc_Flaw = (numreal-nummod)/sigma;
        fRc_Bst2 = (numreal-nummod2)/fStdBst;
        fRc_Bst = (numreal-nummod)/fStdBst;

        % Set line for learning period
        yy = get(gca,'ylim');
        plot([max(time_as) max(time_as)],[0 yy(2)],'k-.')
        string1=['p1 = ' num2str(pval1) '; c1 = ' num2str(cval1) '; k1 = ' num2str(kval1) ];
        string2=['p2 = ' num2str(pval2) '; c2 = ' num2str(cval2) '; k2 = ' num2str(kval2) ];
        string3=['pmed1 = ' num2str(pmed1) '+-' num2str(mStdL(1,1)) '; cmed1 = ' num2str(cmed1) '+-' num2str(mStdL(1,3)) '; kmed1 = ' num2str(kmed1) '+-' num2str(mStdL(1,5))];
        string4=['pmed2 = ' num2str(pmed2) '+-' num2str(mStdL(1,2)) '; cmed2 = ' num2str(cmed2) '+-' num2str(mStdL(1,4)) '; kmed2 = ' num2str(kmed2) '+-' num2str(mStdL(1,6))];
        text(max(time_asf)*0.1,yy(2)*0.9,string1,'FontSize',8);
        text(max(time_asf)*0.1,yy(2)*0.85,string2,'FontSize',8);
        text(max(time_asf)*0.1,yy(2)*0.8,string3,'FontSize',8);
        text(max(time_asf)*0.1,yy(2)*0.75,string4,'FontSize',8);
        string=['\sigma(Bst) = ' num2str(fStdBst) ' Rc(Med) = ' num2str(fRc_Bst) ];%' Rc(Obfit) = ' num2str(fRc_Bst2)];
        text(max(time_asf)*0.1,yy(2)*0.1,string,'FontSize',8);
        %     sInfoStr = ['Omori parameters fitting the original data: ' string1 string2 ...
        %         'Median values from bootstrapping: ' string3 string4];
        %     msgbox(sInfoStr,'Omori law parameters');
        sAIC = ['AIC = ' num2str(fAIC)];
        text(max(time_asf)*0.1,yy(2)*0.05,sAIC,'FontSize',8)
    else
        disp('no result')
    end

author: CelsoReyes


function SM = StreamMaker()
% StreamMaker (Factory for fie/input/output Streams)
%
% A StreamMaker can make Stream PseudoObjects based on either
%  an "fid" or "filename" (and extra arguments for `fopen`).
% The StreamMaker also contains a method `isStream` to validate whether
% the value passed is a valid stream specifier.
%
% Usage
%
%  SM = StreamMaker;
%
%    Stream = SM.make(fid)
%    Stream = SM.make(filename, ...)
%
% This returns a PseudoObject Stream with the following properties:
%   - name: (file) name of the stream
%   - fid: handle (fid) of the stream
%
% and methods:
%   - print: prints to the stream, i.e. fprintf
%   - close: closes the stream, i.e. fclose
%
% It may also contain a field to automatically close the Stream when it goes
% out of scope.
%
    SM = PseudoObject('StreamMaker', ...
                      'isStream',  @isStream, ...
                      'make', @constructStream);
end

function PseudoObj = PseudoObject(T, varargin)
% construct a Pseudo-Object with type T (no other fields yet)
    PseudoObj = struct('Type', T, varargin{:});
end

function bool = isStream(value)
    bool = ischar(value) || ismember(value, [1,2,fopen('all')]);
    %TODO: allow others kinds of streams
    %     Stream -> clipboard (write on close)
    %     Stream -> string variable
    % e.g. a quick-and-dirty way would be to write the file to `tempname`
    % putting a flag to read that file back upon completion.
end

function Stream = constructStream(streamSpecifier, varargin)
    % this is the actual constructor of a stream
    if ~isStream(streamSpecifier)
        error('StreamMaker:NotAStream', 'Invalid stream specifier "%s"', ...
              streamSpecifier);
    end

    Stream = PseudoObject('Stream');
    closeAfterUse = false;
    if ischar(streamSpecifier)
        Stream.name = streamSpecifier;
        Stream.fid = fopen(Stream.name, varargin{:});
        closeAfterUse = true;
    elseif isnumeric(streamSpecifier)
        Stream.fid  = streamSpecifier;
        Stream.name = fopen(Stream.fid);
    end

    if Stream.fid == -1
        error('Stream:InvalidStream', ...
              'Unable to create stream "%s"!', streamSpecifier);
    end

    Stream.print = @(varargin) fprintf(Stream.fid, varargin{:});
    Stream.close = @() fclose(Stream.fid);
    if closeAfterUse
        Stream.closeAfterUse = onCleanup(Stream.close);
    end
end


function A = MetricLearningAutotuneKnn(metric_learn_alg, y, X, params);
% function A = MetricLearningAutotuneKnn(metric_learn_alg, y, X, params);
%
% metric_learn_alg:
% Runs information-theoretic metric learning over various parameters of
% gamma, choosing that with the highest accuracy.
%
% Returns: Mahalanobis matrix A for learned distance metric

if (~exist('params')),
    params = struct();
end
params = SetDefaultParams(params);

% regularize to the identity matrix
A0 = eye(size(X, 2));

% define gamma values for slack variables
gammas = 10.^(-4:4);

accs = zeros(length(gammas), 1);
for (i=1:length(gammas)),
    disp(sprintf('\tTuning burg kernel learning: gamma = %f', gammas(i)));
    params.gamma = gammas(i);
    accs(i) = CrossValidateKNN(y, X, @(y,X) MetricLearning(metric_learn_alg, y, X, A0, params), 2, params.k);
end

[v,i] = max(accs);
gamma = gammas(i);
disp(sprintf('\tOptimal gamma value: %f', gamma));
params.gamma = gamma;
A = MetricLearning(metric_learn_alg, y, X, A0, params);


text:
% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

mdl_puma560

T = p560.fkine(qn);
q = qn;

while true
    T=joy2tr(T, 'tool');
    trprint(T, 'fmt', '%.1f')
    q = p560.ikine6s(T, q);
    p560.plot(q)
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function visualizePts3d(pts3d, R, T, figTitle)
%% Visualize 3d reconstruction results with camera & projector extrinsics
% See also: showExtrinsics

%% License
% ACADEMIC OR NON-PROFIT ORGANIZATION NONCOMMERCIAL RESEARCH USE ONLY
% Copyright (c) 2018 Bingyao Huang
% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:

% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.

% If you publish results obtained using this software, please cite our paper.

% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
% SOFTWARE.

%%

% projector origin in world space
% camOrg = [0,0,0];
Rt = R';
prjOrg = (Rt*(-T))';

% draw
figure(Name=[figTitle ' 3d point cloud'], Units='normalized', Position=[0.3 0.3 0.2 0.3]);
hold on;
title([figTitle ' 3d point cloud'] );

% plot camera and projector models
cameraSize = 20 ;
plotCamera('Size', cameraSize, 'Color', 'b', 'Label', 'Camera', 'Opacity', 0);
grid on

% matlab buildtin function takes rotated OpenCV's R
plotCamera('Location', prjOrg, 'Orientation', Rt', 'Size', cameraSize, ...
    'Color', 'r', 'Label', 'Projector', 'Opacity', 0);

% Label the axes
xlabel('X', 'fontsize', 30);
ylabel('Y', 'fontsize', 30);
zlabel('Z', 'fontsize', 30);

% plot the 3d points in green
if(~isempty(pts3d))
    h = scatter3(pts3d(:,1),pts3d(:,2),pts3d(:,3), 10, 'go', 'filled',...
        'MarkerEdgeColor', 'k');
end

hold off;
daspect([1 1 1]);
view(3);
% axis vis3d tight;
h.Clipping = 'off'; % disable clipping

rotate3d on
end

author: BingyaoHuang

repo: single-shot-pro-cam-calib


text:
clear
close all

% for linux
mex WarpMesh.cpp -lGLU -lOSMesa -lGL
% for mac
%mex WarpMesh.cpp -lGLU -lOSMesa -I/opt/X11/include/ -L/opt/X11/lib/

load demo.mat

%{
P=[570.3422         0  320.0000 0
          0  570.3422  240.0000 0
          0         0    1.0000 0];
%}

XYZcamera = double(XYZcamera);

XYZcamera(:,:,1) = - XYZcamera(:,:,1);
XYZcamera(:,:,2) = XYZcamera(:,:,2);
XYZcamera(:,:,3) = XYZcamera(:,:,3);

[label,depth]=WarpMesh(P,640,480,XYZcamera);

label = label';
label = label(:,end:-1:1);
figure
imagesc(label)

depth = depth';
depth = depth(end:-1:1,end:-1:1);

z_near = 0.3;
z_far_ratio = 1.2;
depth = z_near./(1-double(depth)/2^32);
maxDepth = max(depth(abs(depth) < 100));

cropmask = (depth < z_near) | (depth > z_far_ratio * maxDepth);
depth(cropmask) = 0; %NaN;%z_far_ratio * maxDepth;

figure
imagesc(depth)

author: jianxiongxiao

repo: ProfXkit


text:
function [lp,dlp] = priorGaussMulti(mu,s2,x)

% Multivariate Gaussian hyperparameter prior distribution.
% Compute log-likelihood and its derivative or draw a random sample.
% The prior distribution is parameterized as:
%
%   p(x) = exp(-r2/2) / sqrt(det(2*pi*s2)) where r2(x) = (x-mu)'*inv(s2)*(x-mu),
%
% mu(Dx1) is the mean parameter, s2(Dx1) or s2(DxD) is the variance parameter
% and x(DxN) contains query hyperparameters for prior evaluation.
%
% For more help on design of priors, try "help priorDistributions".
%
% Copyright (c) by Roman Garnett and Hannes Nickisch, 2014-09-12.
%
% See also PRIORDISTRIBUTIONS.M, PRIORGAUSS.M.

if nargin<2, error('mu and s2 parameters need to be provided'), end
if ndims(mu)~=2 || size(mu,2)~=1, error('mu needs to be (Dx1)'), end
D = size(mu,1);
s2_ok = ndims(s2)==2 && all(size(s2)==[D,1] | size(s2)==[D,D]);
if ~s2_ok, error('s2 needs to be (DxD) or (Dx1)'), end
if size(s2,2)==D                                        % full multivariate case
  s = chol(s2)'; lds = sum(log(diag(s)));                 % lds = log(det(s2))/2
else                                                       % diagonal covariance
  s = sqrt(s2);  lds = sum(log(s));
end
if nargin<3                                             % return a random sample
  lp = randn(D,1);                                                 % unit sample
  if size(s,2)==D, lp = s*lp+mu; else lp = s.*lp+mu; end % affine transformation
  return
end
if D==1 && size(x,1)>1                            % both mu/s2 scalar => inflate
  D = size(x,1); mu = mu*ones(D,1); s = s*ones(D,1); lds = D*lds;
end
if ~(ndims(x)==2 && size(x,1)==D), error('x needs to be (Dxn)'), end

oN = ones(1,size(x,2));
if size(s,2)==D
  xs = s\(x-mu*oN);       dlp = -s'\xs;
else
  xs = (x-mu*oN)./(s*oN); dlp = -xs./(s*oN);
end
lp = -sum(xs.^2,1)/2 - D*log(2*pi)/2 - lds;

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [model, varargout] = gmlvq_core(trainSet, trainLab, varargin)
%GMLVQ_core.m - trains the Generalized Matrix LVQ algorithm
%NOTE: minimal requirement version 7.4.0.336 (R2007a)
%
%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  %% Use the wrapper GMLVQ.M to access the functionality in style of %%
%  %% the SOM Toolbox (i.e. with data structs).                       %%
%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  example for usage:
%  trainSet = [1,2,3;4,5,6;7,8,9];
%  trainLab = [1;1;2];
%  GMLVQ_model=GMLVQ_train(trainSet,trainLab); % minimal parameters required
%  estimatedTrainLabels = GMLVQ_classify(trainSet, GMLVQ_model);
%  trainError = mean( trainLab ~= estimatedTrainLabels );
%
% input:
%  trainSet : matrix with training samples in its rows
%  trainLab : vector with the labels of the training set
% optional parameters:
%  PrototypesPerClass: (default=1) the number of prototypes per class used. This could
%  be a number or a vector with the number for each class
%  initialPrototypes : (default=[]) a set of prototypes to start with. If not given initialization near the class means
%  initialMatrix     : the matrix omega to start with. If not given random
%  initialization for rectangular matrices and Unity for squared omega
%  dim               : (default=nb of features for training) the maximum rank or projection dimension
%  regularization    : (default=0) values usually between 0 and 1 treat with care.
%  Regularizes the eigenvalue spectrum of omega'*omega to be more homogeneous
%  testSet           : (default=[]) an optional test set used to compute
%  the test error. The last column is expected to be a label vector
%  comparable        : (default=0) a flag which resets the random generator
%  to produce comparable results if set to 1
%  optimization      : (default=fminlbfgs) indicates which optimization is used: sgd or fminlbfgs
% parameter for the stochastic gradient descent sgd
%  nb_epochs             : (default=100) the number of epochs for sgd
%  learningRatePrototypes: (default=[]) the learning rate for the prototypes.
%  Could be the start and end value used for a sigmoidal spectrum or a vector of length nb_epochs
%  learningRateMatrix    : (default=[]) the learning rate for the matrix.
%  Could be the start and end value used for a sigmoidal spectrum or a vector of length nb_epochs
%  MatrixStart           : (default=1) the epoch to start the matrix training
% parameter for the build-in function fminlbfgs
%  threshstop       : (default=0) the training error for early stopping
%  nb_reiterations  : (default=100) the number of optimization reiterations performed
%  useEarlyStopping : (default=1) use early stopping based on threshstop
%  Display          : (default=iter) the optimization output 'iter' or 'off'
%  GradObj          : (default=on) use the gradient information or not
%  HessUpdate       : (default=lbfgs) the update can be 'lbfgs', 'bfgs' or 'steepdesc'
%  TolFun           : (default=1e-6) the tolerance
%  MaxIter          : (default=2500) the maximal number of iterations
%  MaxFunEvals      : (default=1000000) the maximal number of function evaluations
%  TolX             : (default=1e-10) tolerance
%  DiffMinChange    : (default=1e-10) minimal change
%
% output: the GMLVQ model with prototypes w their labels c_w and the matrix omega
%  optional output:
%  initialization : a struct containing the settings
%  trainError     : error in the training set
%  testError      : error in the training set (only computed if 'testSet' is given)
%  costs          : the output of the cost function
%
% Citation information:
% Petra Schneider, Michael Biehl, Barbara Hammer:
% Adaptive Relevance Matrices in Learning Vector Quantization. Neural Computation 21(12): 3532-3561 (2009)
%
% K. Bunte, P. Schneider, B. Hammer, F.-M. Schleif, T. Villmann and M. Biehl,
% Limited Rank Matrix Learning - Discriminative Dimension Reduction and Visualization,
% Neural Networks, vol. 26, nb. 4, pp. 159-173, 2012.
%
% P. Schneider, K. Bunte, B. Hammer and M. Biehl, Regularization in Matrix Relevance Learning,
% IEEE Transactions on Neural Networks, vol. 21, nb. 5, pp. 831-840, 2010.
%
% Kerstin Bunte (modified based on the code of Marc Strickert http://www.mloss.org/software/view/323/ and Petra Schneider)
% uses the Fast Limited Memory Optimizer fminlbfgs.m written by Dirk-Jan Kroon available at the MATLAB central
% kerstin.bunte@googlemail.com
% Fri Nov 09 14:13:52 CEST 2012
%
% Conditions of GNU General Public License, version 2 and BSD License apply.
% See file 'license-gpl2.txt' and 'BSD_license.txt' enclosed in this package.
% Programs are not for use in critical applications!
%

% Contributed to SOM Toolbox vs2, December 3rd, 2012 by Alexander Schulz
% Copyright (c) Kerstin Bunte
% http://www.cis.hut.fi/projects/somtoolbox/

nout = max(nargout,1)-1;
p = inputParser;   % Create an instance of the class.
p.addRequired('trainSet', @isfloat);
p.addRequired('trainLab', @(x) length(x)==size(trainSet,1) & isnumeric(x));

p.addParamValue('PrototypesPerClass', ones(1,length(unique(trainLab))), @(x)(sum(~(x-floor(x)))/length(x)==1 && (length(x)==length(unique(trainLab)) || length(x)==1)));
p.addParamValue('initialPrototypes',[], @(x)(size(x,2)-1==size(trainSet,2) && isfloat(x)));
p.addParamValue('initialMatrix',[], @(x)(size(x,2)==size(trainSet,2) && isfloat(x)));
p.addParamValue('dim',size(trainSet,2), @(x)(~(x-floor(x)) && x<=size(trainSet,2) && x>0));
p.addParamValue('regularization',0, @(x)(isfloat(x) && x>=0));
p.addOptional('testSet', [], @(x)(size(x,2)-1)==size(trainSet,2) & isfloat(x));
p.addOptional('comparable', 0, @(x)(~(x-floor(x))));
p.addOptional('optimization', 'fminlbfgs', @(x)any(strcmpi(x,{'sgd','fminlbfgs'})));
% parameter for the stochastic gradient descent
p.addOptional('nb_epochs', 100, @(x)(~(x-floor(x))));
p.addParamValue('learningRatePrototypes', [], @(x)(isfloat(x) || isa(x,'function_handle'))); % && (length(x)==2 || length(x)==p.Results.epochs)
p.addParamValue('learningRateMatrix', [], @(x)(isfloat(x)  || isa(x,'function_handle')));
p.addOptional('MatrixStart', 1, @(x)(~(x-floor(x))));
% parameter for the build-in function
p.addOptional('threshstop',0,@(x) isfloat(x));
p.addOptional('nb_reiterations',100,@(x)(~(x-floor(x))));
p.addOptional('useEarlyStopping',1,@(x)(~(x-floor(x))));
p.addOptional('Display', 'off', @(x)any(strcmpi(x,{'iter','off'})));
p.addOptional('GradObj', 'on', @(x)any(strcmpi(x,{'on','off'})));
p.addOptional('HessUpdate', 'lbfgs', @(x)any(strcmpi(x,{'lbfgs', 'bfgs', 'steepdesc'})));
p.addOptional('TolFun',1e-6,@(x) isfloat(x));
p.addOptional('MaxIter', 2500, @(x)(~(x-floor(x))));
p.addOptional('MaxFunEvals', 1000000, @(x)(~(x-floor(x))));
p.addOptional('TolX',1e-10,@(x) isfloat(x));
p.addOptional('DiffMinChange',1e-10,@(x)isfloat(x));
p.CaseSensitive = true;
p.FunctionName = 'GMLVQ';
% Parse and validate all input arguments.
p.parse(trainSet, trainLab, varargin{:});

%%% check if results should be comparable
if p.Results.comparable,
    rng('default');
end
%%% set useful variables
nb_samples = size(trainSet,1);
nb_features = size(trainSet,2);
% labels should be a row vector
if size(trainLab,1)~=nb_samples, trainLab = trainLab';end

classes = unique(trainLab);
nb_classes = length(classes);
dim = p.Results.dim;
MatrixStart = p.Results.MatrixStart;
testSet = p.Results.testSet;
% global regularization;
regularization = p.Results.regularization;
if regularization, disp(['Regularize the eigenvalue spectrum of omega''*omega with ',num2str(regularization)]);end

initialization = rmfield(p.Results, 'trainSet');
initialization.trainSet = [num2str(nb_samples),'x',num2str(nb_features),' matrix'];
initialization = rmfield(initialization, 'trainLab');
initialization.trainLab = ['vector of length ',num2str(length(trainLab))];
if ~isempty(testSet)
    initialization = rmfield(initialization, 'testSet');
    initialization.testSet = [num2str(size(testSet,1)),'x',num2str(size(testSet,2)),' matrix'];
end
switch(p.Results.optimization)
case{'sgd'}
    initialization = rmfield(initialization, 'useEarlyStopping');
    initialization = rmfield(initialization, 'Display');
    initialization = rmfield(initialization, 'GradObj');
    initialization = rmfield(initialization, 'HessUpdate');
    initialization = rmfield(initialization, 'TolFun');
    initialization = rmfield(initialization, 'MaxIter');
    initialization = rmfield(initialization, 'MaxFunEvals');
    initialization = rmfield(initialization, 'TolX');
    initialization = rmfield(initialization, 'DiffMinChange');
    initialization = rmfield(initialization, 'nb_reiterations');
case{'fminlbfgs'}
    disp('The fminlbfgs optimization uses some global variables:');
    disp('threshstop earlystopped useEarlyStopping');
    initialization = rmfield(initialization, 'nb_epochs');
    initialization = rmfield(initialization, 'learningRatePrototypes');
    initialization = rmfield(initialization, 'learningRateMatrix');
    initialization = rmfield(initialization, 'MatrixStart');
    nb_reiterations = p.Results.nb_reiterations;
end
% Display all arguments.
%disp 'Settings for GMLVQ:'
%disp(initialization);

%%% check the number of prototypes per class if one integer is given and turn
%%% it into a vector
nb_ppc = p.Results.PrototypesPerClass;
if length(nb_ppc)~=nb_classes,
    nb_ppc = ones(1,nb_classes)*nb_ppc;
end

%%% initialize the prototypes
if isempty(p.Results.initialPrototypes)
    % initialize near the class centers
    w = zeros(sum(nb_ppc),nb_features);
    c_w = zeros(sum(nb_ppc),1);
    actPos = 1;
    for actClass=1:nb_classes
        nb_prot_c = nb_ppc(actClass);
        classMean = mean(trainSet(trainLab==classes(actClass),:));
        % set the prototypes to the class mean and add a random variation between -0.1 and 0.1
        w(actPos:actPos+nb_prot_c-1,:) = classMean(ones(nb_prot_c,1),:)+(rand(nb_prot_c,nb_features)*2-ones(nb_prot_c,nb_features))/10;
        c_w(actPos:actPos+nb_prot_c-1) = classes(actClass);
        actPos = actPos+nb_prot_c;
    end
else
    % initialize with given w
    w = p.Results.initialPrototypes(:,1:end-1);
    c_w = p.Results.initialPrototypes(:,end);
end
%%% initialize the matrix
if isempty(p.Results.initialMatrix)
    if(p.Results.dim==nb_features)
        omega = eye(nb_features);
    else % initialize with random numbers between -1 and 1
        omega = rand(dim,nb_features)*2-ones(dim,nb_features);
    end
else
    omega = p.Results.initialMatrix;
end
% normalize the matrix
omega = omega / sqrt(sum(diag(omega'*omega)));
model = struct('w',w,'c_w',c_w,'omega',omega);
clear w c_w omega;

% GMLVQ_algorithm = struct('update',@GMLVQ_update,'classify',@GMLVQ_classify,'costfun',@GLVQ_costfun);
switch(p.Results.optimization)
case{'sgd'}
    %%% gradient descent variables
    nb_epochs = p.Results.nb_epochs;
    % compute the vector of nb_epochs learning rates alpha for the prototype learning
    if isa(p.Results.learningRatePrototypes,'function_handle')
        % with a given function specified from the user
        alphas = arrayfun(p.Results.learningRatePrototypes, 1:nb_epochs);
    elseif length(p.Results.learningRatePrototypes)>2
        if length(p.Results.learningRatePrototypes)==nb_epochs
            alphas = p.Results.learningRatePrototypes;
        else
            disp('The learning rate vector for the prototypes does not fit the nb of epochs');
            return;
        end
    else
        % or use an decay with a start and a decay value
        if isempty(p.Results.learningRatePrototypes)
            initialization.learningRatePrototypes = [nb_features/100, nb_features/10000];
        end
        alpha_start = initialization.learningRatePrototypes(1);
        alpha_end = initialization.learningRatePrototypes(2);
        alphas = arrayfun(@(x) alpha_start * (alpha_end/alpha_start)^(x/nb_epochs), 1:nb_epochs);
    %     alphas = arrayfun(@(x) alpha_start / (1+(x-1)*alpha_end), 1:nb_epochs);
    end
    % compute the vector of nb_epochs learning rates epsilon for the Matrix learning
    epsilons = zeros(1,nb_epochs);
    if isa(p.Results.learningRateMatrix,'function_handle')
        % with a given function specified from the user
    % 	epsilons = arrayfun(p.Results.learningRateMatrix, 1:nb_epochs);
        epsilons(MatrixStart:nb_epochs) = arrayfun(p.Results.learningRateMatrix, MatrixStart:nb_epochs);
    elseif length(p.Results.learningRateMatrix)>2
        if length(p.Results.learningRateMatrix)==nb_epochs
            epsilons = p.Results.learningRateMatrix;
        else
            disp('The learning rate vector for the Matrix does not fit the nb of epochs');
            return;
        end
    else
        % or use an decay with a start and a decay value
        if isempty(p.Results.learningRateMatrix)
            initialization.learningRateMatrix = [nb_features/1000, nb_features/100000];
        end
        eps_start = initialization.learningRateMatrix(1);
        eps_end = initialization.learningRateMatrix(2);
    %     epsilons = arrayfun(@(x) eps_start * (eps_end/eps_start)^(x/nb_epochs), 1:nb_epochs);
        epsilons(MatrixStart:nb_epochs) = arrayfun(@(x) eps_start * (eps_end/eps_start)^((x-MatrixStart)/(nb_epochs-MatrixStart)), MatrixStart:nb_epochs);
    end

    %%% initialize requested outputs
    trainError = [];
    costs = [];
    testError = [];
    if nout>=2,
        % train error requested
        trainError = ones(1,nb_epochs+1);
        estimatedLabels = GMLVQ_classify(trainSet, model); % error after initialization
        trainError(1) = sum( trainLab ~= estimatedLabels )/nb_samples;
        if nout>=3,
            % test error requested
            if isempty(testSet)
                testError = [];
                disp('The test error is requested, but no labeled test set given. Omitting the computation.');
            else
                testError = ones(1,nb_epochs+1);
                estimatedLabels = GMLVQ_classify(testSet(:,1:end-1), model); % error after initialization
                testError(1) = sum( testSet(:,end) ~= estimatedLabels )/length(estimatedLabels);
            end
            if nout>=4,
                % costs requested
%                 LabelEqPrototype = trainLab*ones(1,numel(model.c_w)) == (model.c_w*ones(1,nb_samples))';
                disp('The computation of the costs is an expensive operation, do it only if you really need it!');
                costs = ones(1,nb_epochs+1);
                costs(1) = GMLVQ_costfun(trainSet, trainLab, model, regularization);
%                 costs(1) = sum(arrayfun(@(idx) GLVQ_costfun(min(dist(idx,model.c_w == trainLab(idx))),...
%                                                             min(dist(idx,model.c_w ~= trainLab(idx))))-regTerm, 1:size(dist,1)));
            end
        end
    end

    %%% optimize with stochastic gradient descent
    for epoch=1:nb_epochs
        if mod(epoch,100)==0, disp(epoch); end
        % generate order to sweep through the trainingset
        order = randperm(nb_samples);
        % perform one sweep through trainingset
        for i=1:nb_samples
            % select one training sample randomly
            xi = trainSet(order(i),:);
            c_xi = trainLab(order(i));

%             dist = ((xi(ones(size(model.w,1),1),:))-model.w)*model.omega'*(model.omega*((xi(ones(size(model.w,1),1),:))-model.w)');
%             dist = diag(dist);
            dist = sum((bsxfun(@minus, xi, model.w)*model.omega').^2, 2);
            % determine the two winning prototypes
            % nearest prototype with the same class
            [sortDist,sortIdx] = sort(dist);
            count = 1;
            J = sortIdx(count);
            while model.c_w(sortIdx(count)) ~= c_xi,
                count = count+1;
                J = sortIdx(count);
            end
            dJ = sortDist(count);
            count = 1;
            K = sortIdx(count);
            while model.c_w(sortIdx(count)) == c_xi,
                count = count+1;
                K = sortIdx(count);
            end
            dK = sortDist(count);
    %         disp([J,K,dJ,dK]);
            wJ = model.w(J,:);
            wK = model.w(K,:);
            % prototype update
            norm_factor = (dJ + dK)^2;
            DJ = (xi-wJ);
            DK = (xi-wK);

            oo = model.omega'*model.omega;

            dwJ = (2*dK/norm_factor)*2*oo*DJ';
            dwK = (2*dJ/norm_factor)*2*oo*DK';
    %         dwJ = (2*dK/norm_factor)*2*model.omega'*model.omega*DJ';
    %         dwK = (2*dJ/norm_factor)*2*model.omega'*model.omega*DK';
            model.w(J,:) = wJ + alphas(epoch) * dwJ';
            model.w(K,:) = wK - alphas(epoch) * dwK';
            % update matrices
            if epsilons(epoch)>0, % epoch >= MatrixStart
                f1 = (2*dK/norm_factor)*2*(model.omega*DJ')*DJ;
                f2 = (2*dJ/norm_factor)*2*(model.omega*DK')*DK;
                % update omega
                if regularization,
                    f3 = (pinv(model.omega))';
                else
                    f3 = 0;
                end
                model.omega = model.omega-epsilons(epoch) * (f1-f2  - regularization * f3);
                % normalization
                model.omega = model.omega / sqrt(sum(diag(oo)));
            end
        end
        if nout>=2,
            % train error requested
            estimatedLabels = GMLVQ_classify(trainSet, model); % error after epoch
            trainError(epoch+1) = sum( trainLab ~= estimatedLabels )/nb_samples;
            if nout>=3,
                % test error requested
                if ~isempty(testSet)
                    estimatedLabels = GMLVQ_classify(testSet(:,1:end-1), model); % error after initialization
                    testError(epoch+1) = sum( testSet(:,end) ~= estimatedLabels )/length(estimatedLabels);
                end
                if nout>=4,
                    % costs requested
                    costs(epoch+1) = GMLVQ_costfun(trainSet, trainLab, model, regularization);
%                     costs(epoch+1) = sum(arrayfun(@(idx) GLVQ_costfun(min(dist(idx,model.c_w == trainLab(idx))),...
%                                                                       min(dist(idx,model.c_w ~= trainLab(idx))))-regTerm, 1:size(dist,1)));
                end
            end
        end
    end
case{'fminlbfgs'}
    %%% optimization options
    options = struct( ...
      'Display',p.Results.Display, ...
      'GradObj',p.Results.GradObj, ...
      'GradConstr',false, ...
      'GoalsExactAchieve',0, ...
      'TolFun',p.Results.TolFun, ...
      'MaxIter',p.Results.MaxIter, ...
      'MaxFunEvals', p.Results.MaxFunEvals, ...
      'TolX',p.Results.TolX, ...
      'DiffMinChange',p.Results.DiffMinChange, ...
      'OutputFcn','LVQ_progresser', ...
      'HessUpdate',p.Results.HessUpdate ...
    );
    clear('progresser'); % memory therein might need reset
    global threshstop earlystopped useEarlyStopping % for LVQ_progresser.m datval labval n_vec
    useEarlyStopping = p.Results.useEarlyStopping; % use early stopping
    earlystopped = false;
    threshstop = p.Results.threshstop; % stop if classification below this threshold for early stopping
%     prototypeLabel = model.c_w;
    nb_prototypes = numel(model.c_w);
%     LabelEqualsPrototype = trainLab*ones(1,nb_prototypes) == (model.c_w*ones(1,nb_samples))';
    LabelEqualsPrototype = bsxfun(@eq,trainLab,model.c_w');
    earlystopped = false; % don't change, assigned in progresser.m
    clear('progresser'); % memory therein might need reset
    newfval = realmax('single');
    % fminlbfgs optimizer courtesy of Dirk-Jan Kroon:
    % http://www.mathworks.de/matlabcentral/fileexchange/23245
    % early stopping to be implemented in progresser function
    variables = zeros(dim+nb_prototypes,size(trainSet,2));
    variables(1:nb_prototypes,:) = model.w;
    variables(nb_prototypes+1:end,:) = model.omega;
    LRprototypes = 1; % learn prototype locations
    LRrelevances = 0; % don't learn metric
%     [variables,fval] = fminlbfgs(@GMLVQ_optfun,variables,options);
    [variables,fval] = fminlbfgs(@(variables) GMLVQ_optfun(variables,trainSet,LabelEqualsPrototype,LRrelevances,LRprototypes,model.c_w,regularization),variables,options);
    if not(isempty(fval))
        newfval = fval;
    end
    LRprototypes = 0; % don't learn prototype locations
    LRrelevances = 1; % learn metric
%     [variables,fval] = fminlbfgs(@GMLVQ_optfun,variables,options);
    [variables,fval] = fminlbfgs(@(variables) GMLVQ_optfun(variables,trainSet,LabelEqualsPrototype,LRrelevances,LRprototypes,model.c_w,regularization),variables,options);

    if not(isempty(fval))
        newfval = fval;
    end
    if not(isempty(fval))
      clear('progresser'); % memory therein might need reset
      LRprototypes = 1;
      for i = 1:nb_reiterations  % depending on data, re-iterations might further improve
%         [variables,fval] = fminlbfgs(@GMLVQ_optfun,variables,options);
        [variables,fval] = fminlbfgs(@(variables) GMLVQ_optfun(variables,trainSet,LabelEqualsPrototype,LRrelevances,LRprototypes,model.c_w,regularization),variables,options);
        if not(isempty(fval))
          if abs(fval - newfval) < 1e-3
            newfval = fval;
            break
          end
          newfval = fval;
        end
        if isempty(fval) || earlystopped
          break
        end
      end
    end
    model.w = variables(1:nb_prototypes,:);
    model.omega = variables(nb_prototypes+1:end,:);
    model.omega = model.omega / sqrt(sum(diag(model.omega'*model.omega)));
    if nout>=2,
        % train error requested
        estimatedLabels = GMLVQ_classify(trainSet, model); % error after initialization
        trainError = mean( trainLab ~= estimatedLabels );
        if nout>=3,
            % test error requested
            if isempty(testSet)
                testError = [];
                disp('The test error is requested, but no labeled test set given. Omitting the computation.');
            else
                estimatedLabels = GMLVQ_classify(testSet(:,1:end-1), model); % error after initialization
                testError = mean( testSet(:,end) ~= estimatedLabels );
            end
            if nout>=4,
                % costs requested
                costs = GMLVQ_costfun(trainSet, trainLab, model, regularization);
%                 dist = computeDistance(trainSet, model.w, model);
%                 if regularization,
%                     regTerm = regularization * log(det(model.omega*model.omega'));
%                 else
%                     regTerm = 0;
%                 end
%                 costs = sum(arrayfun(@(idx) GLVQ_costfun(min(dist(idx,model.c_w == trainLab(idx))),...
%                                                          min(dist(idx,model.c_w ~= trainLab(idx)))), 1:size(dist,1)))-regTerm;
            end
        end
    end
end
%%% output of the training
varargout = cell(nout);
for k=1:nout
	switch(k)
		case(1)
			varargout(k) = {initialization};
		case(2)
			varargout(k) = {trainError};
		case(3)
			varargout(k) = {testError};
		case(4)
            varargout(k) = {costs};
	end
end

function cost = GMLVQ_costfun(trainSet, trainLab, model, regularization)
%GMLVQ_costfun.m - computes the costs for a given training set and GMLVQ
%model with or without regularization
%  example for usage:
%  trainSet = [1,2,3;4,5,6;7,8,9];
%  trainLab = [1;1;2];
%  GMLVQ_model=GMLVQ_train(trainSet,trainLab); % minimal parameters required
%  costs = GMLVQ_costfun(trainSet, trainLab, GMLVQ_model, 0);
%
% input:
%  trainSet : matrix with training samples in its rows
%  trainLab : a vector of training labels
%  model    : GMLVQ model with prototypes w their labels c_w and the matrix omega
%  regularization: the factor>=0 for the regularization
%
% output    : cost function value
%
% Kerstin Bunte (based on the code from Marc Strickert)
% kerstin.bunte@googlemail.com
% Mon Nov 05 09:05:52 CEST 2012
%
% Conditions of GNU General Public License, version 2 apply.
% See file 'license-gpl2.txt' enclosed in this package.
% Programs are not for use in critical applications!
%
nb_samples = length(trainLab);
% labels should be a row vector
if size(trainLab,1)~=nb_samples, trainLab = trainLab';end

% LabelEqPrototype = trainLab*ones(1,numel(model.c_w)) == (model.c_w*ones(1,nb_samples))';
LabelEqPrototype = bsxfun(@eq,trainLab,model.c_w');
dists = computeDistance(trainSet, model.w, model);
if regularization,
    regTerm = regularization * log(det(model.omega*model.omega'));
% if strcmp(p.Results.optimization,'sgd')
else
    regTerm = 0;
end
Dwrong = dists;
Dwrong(LabelEqPrototype) = realmax(class(Dwrong));   % set correct labels impossible
distwrong = min(Dwrong.'); % closest wrong
clear Dwrong;

Dcorrect = dists;
Dcorrect(~LabelEqPrototype) = realmax(class(Dcorrect)); % set wrong labels impossible
distcorrect = min(Dcorrect.'); % closest correct
clear Dcorrect;
clear dists;
distcorrectpluswrong = distcorrect + distwrong;
distcorrectminuswrong = distcorrect - distwrong;
mu = distcorrectminuswrong ./ distcorrectpluswrong;
if regularization,
    regTerm = regularization * log(det(model.omega*model.omega'));
else
    regTerm = 0;
end
cost = sum(mu)-regTerm;

function [f G]  = GMLVQ_optfun(variables,training_data,LabelEqualsPrototype,LRrelevances,LRprototypes,prototypeLabel,regularization)
% [f G] = GMLVQ_optfun(variables)
% function to be optimzed by matrix relevance learning vector quantization
% variables = [prototype matrix;omega matrix]
% global variables are
%   training_data        : data vectors as row vectors, i.e. attributes in columns
%   LabelEqualsPrototype : binary matrix indicating coocurrences of
%   training labels and prototype labels
%   prototypeLabel       : label vector for the prototypes
%   regularization       : the regularization parameter
%   LRrelevances         : learning rate for the relevance matrix
%   LRprototypes         : learning rate for the prototypes
%
% Kerstin Bunte (modified based on the code of Marc Strickert http://www.mloss.org/software/view/323/)
% kerstin.bunte@googlemail.com
% Fri Nov 09 14:13:52 CEST 2012
%
% Conditions of GNU General Public License, version 2 apply.
% See file 'license-gpl2.txt' enclosed in this package.
% Programs are not for use in critical applications!
%
if isempty(LRprototypes) % values between 1e-2,1e-3,... 1e-8 seem pragmatic
    LRprototypes = 1; % no relevance learning by default
end
if isempty(LRrelevances) % values between 1e-2,1e-3,... 1e-8 seem pragmatic
    LRrelevances = 0; % no relevance learning by default
end
[n_data, n_dim] = size(training_data);
nb_prototypes =  numel(prototypeLabel);
omegaT = variables(nb_prototypes+1:end,:)';
n_vec = size(variables,1) - nb_prototypes;

dists = squaredEuclidean(training_data*omegaT, variables(1:nb_prototypes,:)*omegaT);

Dwrong = dists;
Dwrong(LabelEqualsPrototype) = realmax(class(Dwrong));   % set correct labels impossible
[distwrong pidxwrong] = min(Dwrong.'); % closest wrong
clear Dwrong;

Dcorrect = dists;
Dcorrect(~LabelEqualsPrototype) = realmax(class(Dcorrect)); % set wrong labels impossible
[distcorrect pidxcorrect] = min(Dcorrect.'); % closest correct
clear Dcorrect;

distcorrectpluswrong = distcorrect + distwrong;
distcorrectminuswrong = distcorrect - distwrong;
mu = distcorrectminuswrong ./ distcorrectpluswrong;
% callitq = 1./(1 + exp(-squashsigmoid * mu)); % apply sigmoidal

if regularization,
    regTerm = regularization * log(det(omegaT'*omegaT));
else
    regTerm = 0;
end
f = sum(mu)-regTerm;
% f = mean(callitq);

if nargout > 1  % gradient needed not just function eval
    G = zeros(size(variables)); % initially no gradient
    %       callitq = squashsigmoid * callitq .* (1-callitq); % derivative of sigmoid
    %       distcorrectpluswrong = 2 * callitq ./ distcorrectpluswrong.^2; % degeneration?
    distcorrectpluswrong = 4 ./ distcorrectpluswrong.^2; % norm_factor for derivative for every data sample
    if LRrelevances > 0
        Gw = zeros(n_vec,n_dim);
    end
    for k=1:nb_prototypes%(n_vec+1):size(lambda,1) % update all prototypes
        idxc = (k == pidxcorrect);  % Js: idxs where actual prototype is nearest correct
        idxw = (k == pidxwrong);    % Ks: idxs where actual prototype is nearest wrong

        dcd =  distcorrect(idxw) .* distcorrectpluswrong(idxw);
        dwd =    distwrong(idxc) .* distcorrectpluswrong(idxc);
        if LRrelevances > 0
            % part of derivative of distance
            difc = bsxfun(@minus,training_data(idxc,:),variables(k,:)); % DJs
            difw = bsxfun(@minus,training_data(idxw,:),variables(k,:)); % DKs
            % update omega
            Gw = Gw - (bsxfun(@times,difw,dcd.') * omegaT).' * difw + ...
                      (bsxfun(@times,difc,dwd.') * omegaT).' * difc;
            if LRprototypes > 0
                G(k,:) = dcd * difw - dwd * difc;
            end
        else
            if LRprototypes > 0
                G(k,:) = dcd * training_data(idxw,:) - dwd * training_data(idxc,:) + (sum(dwd)-sum(dcd)) * variables(k,:);
            end
        end
    end
if regularization,
    f3 = (pinv(omegaT'))';
else
    f3 = 0;
end
    % some rescalings needed
    if LRrelevances > 0
        G(nb_prototypes+1:nb_prototypes+n_vec,:) = 2/n_data * LRrelevances * Gw - regularization*f3;
    end
    if LRprototypes > 0
        G(1:nb_prototypes,:) = 1./n_data * LRprototypes * G(1:nb_prototypes,:) * omegaT * omegaT.';
    end
    G = G .* (1 + .0001 * (rand(size(G))-.5)); % help break symmetries
end
% if 0,
% w = variables(1:nb_prototypes,:);
% dJs = zeros(1,nb_samples);
% dKs = zeros(1,nb_samples);
% Js = zeros(1,nb_samples);
% Ks = zeros(1,nb_samples);
% norm_factors = zeros(1,nb_samples);
% DJs = zeros(nb_samples,size(training_data,2));
% DKs = zeros(nb_samples,size(training_data,2));
% for i=1:nb_samples
%     % select one training sample randomly
%     xi = training_data(i,:);
%     c_xi = trainLab(i);
%
%     dist = ((xi(ones(size(w,1),1),:))-w)*omegaT*(omegaT'*((xi(ones(size(w,1),1),:))-w)');
%     dist = diag(dist);
%     % determine the two winning prototypes
%     % nearest prototype with the same class
%     [sortDist,sortIdx] = sort(dist);
%     count = 1;
%     J = sortIdx(count);
%     while prototypeLabel(sortIdx(count)) ~= c_xi,
%         count = count+1;
%         J = sortIdx(count);
%     end
%     dJ = sortDist(count);
%     dJs(i) = dJ;
%     Js(i) = J;
%     count = 1;
%     K = sortIdx(count);
%     while prototypeLabel(sortIdx(count)) == c_xi,
%         count = count+1;
%         K = sortIdx(count);
%     end
%     dK = sortDist(count);
%     dKs(i) = dK;
%     Ks(i) = K;
%
%     wJ = w(J,:);
%     wK = w(K,:);
%     % prototype update
%     norm_factors(i) = 4/((dJ + dK)^2);
%     DJ = (xi-wJ);
%     DK = (xi-wK);
%     DJs(i,:) = DJ;
%     DKs(i,:) = DK;
% end
% end

function D = squaredEuclidean(A, B)
% computes the sqared Euclidean distance
%
%   D = squaredEuclidean(X) returns the squared Euclidean distance matrix of data in rows of X
%   D = squaredEuclidean(X, Y) returns the distance matrix with all distances between the points in X and Y.
%
if nargin == 1 % means that one matrix
    D = bsxfun(@plus, sumsquared(A,2), bsxfun(@minus, sumsquared(A,2).', 2*A*A.')); % 2*(Y*Y.')
else
    D = bsxfun(@plus, sumsquared(A,2), bsxfun(@minus, sumsquared(B,2).', 2*A*B.'));
end
D = max(D,0);

function sq  = sumsquared(x,dim)
% sq  = sumsquared(x,dim)
% sum of all squared element of matrix x along dimension dim

persistent isoctave

if isempty(isoctave)
  isoctave = exist('OCTAVE_VERSION','builtin');
end

if nargin == 1
  dim = 1;
end

if isoctave
  sq = sumsq(x,dim);
else
  sq = sum(x.^2, dim);
end

author: ilarinieminen

repo: SOM-Toolbox


text:
classdef SSC_NN_OMP < handle
    % Implements sparse subspace clustering algorithm using OMP algorithm
    properties
        Quiet = false
    end

    properties(SetAccess=private)
        % Dimension of signal space
        N
        % Number of signals in the space
        S
        % A matrix of size NxS where each column is one signal vector
        Data
        % data vectors after normalization
        NormalizedData
        % The sparsity level or the largest dimension of the sparse subspaces
        K
        % The expected number of subspaces
        NumSubspaces
        % Labels obtained through clustering
        Labels
        % Representation matrix (each column is a representation vector for one data vector)
        Representation
        % affinities identified during sparse coding
        Affinity
        % Adjacency matrix
        Adjacency
        % The spectral clusterer used in the algorithm
        Clusterer
    end

    methods
        function self = SSC_NN_OMP(X, K, NumSubspaces)
            % Constructor
            self.Data = X;
            self.K = K;
            if nargin < 3
                NumSubspaces = -1;
            end
            self.NumSubspaces = NumSubspaces;
            [n, s] = size(X);
            self.N = n;
            self.S = s;
            self.Labels = ones(s, 1);
            % Each data vector is represented using other data vectors in same space
            self.Representation = zeros(s, s);
            self.Affinity = zeros(s, s);
        end

        function result = solve(self)
            % prepare sparse representations
            self.recover_coefficients();
            self.build_adjacency();

            % conduct spectral clustering
            options.num_clusters = self.NumSubspaces;
            result = spx.cluster.spectral.simple.normalized_symmetric(self.Adjacency, options);
            cluster_labels = result.labels;

            % We are disabling our version of spectral clustering for now.
            % clusterer = spx.cluster.spectral.Clustering(self.Adjacency);
            % keep reference for debugging purposes.
            % self.Clusterer = clusterer;
            % clusterer.NumClusters = self.NumSubspaces;
            %cluster_labels = clusterer.cluster_random_walk();

            self.Labels = cluster_labels;
            % Return the labels
            result.Labels = self.Labels;
            result.Z = self.Representation;
            result.W = self.Adjacency;

        end

    end

    methods(Access=private)

        function recover_coefficients(self)
            % Computes sparse representations of the data vectors
            data = self.Data;
            self.NormalizedData = spx.norm.normalize_l2(data);
            % Number of data vectors
            ns = self.S;
            % iterate over signals
            all_cols = 1:ns;
            for s=all_cols
                if ~self.Quiet
                    fprintf('.');
                    if mod(s, 50) == 0
                        fprintf('\n');
                    end
                end
                self.nn_omp(s);
            end
            if ~self.Quiet
                fprintf('\n');
            end
        end

        function nn_omp(self, s)
            % Apply nearest neighbor OMP on s-th vector.
            % data matrix
            X = self.NormalizedData;
            % Current vector
            x =X(:, s);
            % number of vectors
            ns = self.S;
            % Subspace dimension
            nk = self.K;
            % current residual
            r = x;
            % norm of current residual
            res_norm = norm(r);
            % ambient dimension
            nd = self.N;
            % coefficients
            z  = zeros(ns, 1);
            % affinities
            c = zeros(ns, 1);
            % selected indices
            omega = [];
            % max number of iterations
            maxIter = nk;
            MaxResNorm = 1/8;
            % maximum number of nearest neighbors
            nn = min (max(2*(nk-1), 2), round(ns/2) );
            for iter=1:maxIter
                % identify normalized inner product
                rr = r / res_norm;
                % Compute inner products with normalized residual
                inner_products = X' * rr;
                inner_products(s) = 0;
                inner_products(omega) = 0;
                inner_products = abs(inner_products);
                % Find the highest inner product
                [sorted_inner_products, indices] = sort(inner_products, 'descend');
                index = indices(1);                % Add this index to support
                omega = [omega, index];
                angles  = rad2deg(acos(sorted_inner_products(1:nn)'));
                if false
                    fprintf('%.0f ', angles);
                    fprintf('\n');
                end
                min_angle = angles(1);
                % allow for some degree leeway
                max_angle = min(min_angle + 6, 1.5 * min_angle);
                max_index = find(angles > max_angle, 1);
                if isempty(max_index)
                    % all vectors are too close by
                    max_index = nn;
                end
                % number of neighbors
                nnn = max(max_index - 1, 1);
                %fprintf('%d:%d ', nn, nnn);
                c(indices(1:nnn)) = inner_products(indices(1:nnn));
                % Solve least squares problem
                subdict = X(:, omega);
                tmp = linsolve(subdict, x);
                % Updated solution
                z(omega) = tmp;
                % Let us update the residual.
                r = x - subdict * tmp;
                res_norm = norm(r);
                if res_norm < MaxResNorm
                    break;
                end
            end
            self.Representation(:, s) = z;
            self.Affinity(:, s) = c;
        end

        function detect_outliers(self)
            % Identifies outliers in the representations
        end

        function build_adjacency(self)
            C = abs(self.Affinity);
            % Normalize the matrix by column wise maximums
            C = spx.norm.normalize_linf(C);
            % Make it symmetric
            C = C + C';
            % Keep it
            self.Adjacency = C;
        end

    end

end

author: indigits

repo: sparse-plex


clear all
% ===== make testdata =====
% ti=0;
% tf=10;
% dt=.01;
% t=ti:dt:tf; % second
%
% a1=1;
% f1=5;
% y1=a1*cos(2*pi*f1*t)*rand(length(t));
% a2=5;
% f2=1;
% y2=a2*sin(2*pi*f2*t);
%
% testdata=y1+y2+10;

filename = 'cmy(finish).wav';
[data,Fs]=audioread(filename);

% wavplay(y,Fs)
testdata=data';
dt=1/Fs;
t=0:dt:length(testdata)/Fs-dt;

% ===== decomposition in time domain =====
st=100; % sifting time; 100-300
[IMF]=EMD(testdata,st);
% save('IMF.mat','IMF')
% load IMF.mat

%  for i=1:10;
%      plot(IMF(i,:));
%      title(['IMF' num2str(i)])
%      pause
%      close all
%  end

% set(gca,'xtick',[0:100:1000],'xticklabel',0:10);
% set(gca,'ytick',get(gca,'ytick'),'yticklabel',{'DC term','IMF5','IMF4','IMF3','IMF2','IMF1','test data'});
% xlabel('Time (second)');
%  pause
% ===== computation of instantaneous frequency and amplitude =====
[ipifs,ifs,amp]=imf2ipfa(IMF(1:end-1,:),dt); % the DC term should not be include

% ===== Hilbert spectrum =====
% x=round([t;t;t;t;t]'*100+1); % xlabel, time
[L,W]=size(ifs);

x=NaN(L,W);
for i=1:W;
    x(:,i)=t';
end
y=ifs;

x=x(:);
y=y(:);
z=amp(:);

% === x-axis (time) ===
dx=0.01;% time resolution 0.1 second
xi=0;
xf=t(end);

T=[xi:dx:0.5];

% === y-axis (frequency) ===
dy=10;% frequency resolution 100 Hz
yi=0;
yf=1500;
F=yi:dy:yf;

z(y>yf)=[];
x(y>yf)=[];
y(y>yf)=[]; % romove the data with frequency > 5000Hz

x=round(x/dx)+1;
y=round(y/dy)+1; % ylabel, frequency

x(z<0)=[]; % remove the bad data (negative amplitude)
y(z<0)=[];
z(z<0)=[];

% pause

% z(y>90)=[];
% x(y>90)=[];
% y(y>90)=[]; % romove the data with frequency > 9Hz

m=zeros(length(F),length(T));
c=m;
for i=1:length(z);
    m(y(i),x(i))=m(y(i),x(i))+z(i);
    c(y(i),x(i))=c(y(i),x(i))+1;
end

hsp=m./c;

close all
figure('position',[300 100 600 600])
axes('position',[.15 .64 .7 .3])
plot(t,testdata,'-k');
title('(a) Data');
set(gca,'fontsize',16)
grid on
axis tight

axes('position',[.15 .12 .7 .4])
pcolor(T,F,hsp);
shading flat
%[C,H]=contourf(T,F,hsp);
%set(H,'levellist',[0:.0001:.3])
caxis([0 .02])
cb=colorbar('position',[.88 .12 .03 .4]);
set(get(cb,'title'),'string','Amplitude')
xlabel('Time (second)');
ylabel('Frequency (Hz)');
ylim([0 1500])
title('(b) Hilbert spectrum');
set(gca,'fontsize',16)
set(gcf,'paperpositionmode','auto')
print -dpng -r300 Fig01


text:
function class_calc = knnclass2(neighbors_class,neighbors_distance,num_class,K)

% find class on the basis of K neighbors
for g = 1:num_class
    freq(g) = length(find(neighbors_class == g));
end

unique_class = find(freq == max(freq));

if neighbors_distance(1)<1e-5
    class_calc=neighbors_class(1);
elseif length(unique_class) == 1
    [M,class_calc] = max(freq);
else
    mean_dist = ones(num_class,1).*max(neighbors_distance);
    for g = 1:length(unique_class)
        in = find(neighbors_class == unique_class(g));
        mean_dist(unique_class(g)) = mean(neighbors_distance(in));
    end
    [m,class_calc] = min(mean_dist);
end

author: kmansouri

repo: OPERA

save_path: github-repos/MATLAB/kmansouri-OPERA


text:
%% 2D Laplace Equation
% This program solves the following PDE:
%
% $$\frac{\partial^2{T}}{\partial{x^2}}+\frac{\partial^2{T}}{\partial{y^2}}=0$$

clear
clc
close all

%% Discretization of space
% Length and heigth of our 2D domain:
L = 1; % length in x direction
H = 1; % height in y direction

% our target is to evaluate the step h size, thus
h = 0.05; % we asume and equally spaced XY grid
hx = h;
hy = h;

% number of points
n = round(L/hx);
m = round(H/hy);

% initaliza Matriz T
T=zeros(n+1,m+1);

%% Boundary Contidions
% Our domain requires 4 boudary conditions in total, can be only nuemann,
% or only dirichlet or a combination of both.

%% Nuemann Boundary Conditions (Prescribed values at the borders)
% Bottom (last row of the matrix)
T(1,:)=0;
% Top (first row of the matrix)
%T(1,:)=1;
for i=0:n
    %using function of x to prescribe the values.
    T(m+1,i+1)=1-cos(pi/(n)*i);
end
% Left
%T(:,1)=0;
% Right
%T(:,n+1)=2;

%T

% Vector of Prescribed Bondary Values
% Bottom
fb=(T(1,2:n)).' ; %transpose this array
% Top
ft=T(m+1,2:n).' ; %transpose this array
% Left
fl=T(2:n,1);
% Right
fr=T(2:n,n+1);

% Computing vector of known values
fs=zeros(n-1,1);
fs(1)=fl(1);
fs(n-1)=fr(1);
w=fb+fs;
for i=2:n-2
    fs(1)=fl(i);
    fs(n-1)=fr(i);
    w=cat(1,w,fs);
end
fs=zeros(n-1,1);
fs(1)=fl(n-1);
fs(n-1)=fr(n-1);
y=ft+fs;
w=cat(1,w,y);

%% Dirichlect Boundary Conditions (flux at the borders of the domain)
% This conditions requires de mofification of the diagonal values of our
% solution matrix:

% Bottom (last row of the matrix)
%sb=ones(n-1,1); % to insulate this boundary
sb=zeros(n-1,1); %to suppres this condition
% Top (first row of the matrix)
%st=ones(n-1,1); % to insulate this boundary
st=zeros(n-1,1); %to suppres this condition
% Left
sl=zeros(m-1,1);
sl(1)=1; % to insulate this boundary
% Right
sr=zeros(m-1,1);
sr(m-1)=1; % to insulate this boundary

% Computing the diagonal matrix to substract
s=sl+sr; %sum condition in the right and left
stt=st+s; %sum condition of the top and sides
sbb=sb+s; %sum condition of the bottom and sides
u=cat(1,s,s);
for i=4:n-2
    u=cat(1,u,s);
end
z=cat(1,stt,u);
z=cat(1,z,sbb);
Z=spdiags(z,0,length(z),length(z)); %transforms the vector in a diagonal mat.
%spy(Z)

%% System Solution
% Our solution has the shape:
%
% $$ A_{ij}u_j=f_i  $$

% Formulating A:
I=eye(m-1);
e=ones(m-1,1);
K=spdiags([e,-4*e,e],[-1,0,1],m-1,m-1);
S=spdiags([e,e],[-1,1],m-1,m-1);
A=(kron(I,K)+kron(S,I));
A=(A+Z); %(for flux conditions)
u=inv(A)*(-w);
spy(A)

%% Re-arranging the data in vector u

for j=1:n-1
    for i=1:n-1
        t(j,i)=u(i+(n-1)*(j-1));
    end
end
for j=1:n-1
    for i=1:n-1
        T(i+1,j+1)=t(i,j);
    end
end

% Make up for the plot
if sum(sb)>=1
    T(1,:)=T(2,:);
end
if sum(st)>=1
    T(n+1,:)=T(n,:);
end
if sum(sl)>=1
    T(:,1)=T(:,2);
end
if sum(sr)>=1
    T(:,n+1)=T(:,n);
end
%T
%surf(T)
contourf(T)
colormap hot
colorbar('location','southoutside')

%% Computing the Error of our aproximation
% % The exact solution of our problem is given by:
% %
% % $$ T(x,y)=y-\frac{sinh(\pi y))}{sinh(\pi)}cos(\pi x) $$
% for j=1:m+1
%     for i=1:n+1
%         T2(j,i)=j/(m+1)-sinh(pi/(m+1)*j)/sinh(pi)*cos(pi/(n+1)*i);
%     end
% end
% Error=T2-T;
% %T2
% h
% NormError=norm(Error)
% %surf(T2)
% contourf(T2)
% colormap hot
% colorbar('location','southoutside')

repo: Aero-matlab


text:
% Sign is the enumeration class for negative and positive signs.
classdef Sign < Enumerated
	enumeration
		n('negative')
		p('positive')
	end

	methods (Static)
		function elems = elems(ind)
			elems = [Sign.n, Sign.p];
			if nargin > 0  % ind
				elems = elems(ind);
			end
		end

		function count = count()
			count = length(Sign.elems);
		end
	end

	methods
		function g = alter(this)
			% This function simply returned Sign.p if this == Sign.n and Sign.n
			% if this == Sign.p originally, but it is modified to handle
			% arguments given as arrays.
			chkarg(istypeof(this, 'Sign'), '"this" should have Sign as elements.');
			if isempty(this)  % alter([]) makes no sense, but alter(Sign.empty()) does
				g = [];
			else
				g = this.elems(GT.count + 1 - int(this));
			end
		end
	end
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
% The COBRAToolbox: testMets2str.m
%
% Purpose:
%     - testMets2str tests the functionality of mets2str in the
%     rBioNet extension
%
% Authors:
%     - Stefania Magnusdottir April 2017
%

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testMets2str'));
cd(fileDir);

% load rBioNet reaction database
load([fileDir, filesep 'rxn.mat'])

% concatenate cell array of strings into one string
outstr = mets2str(rxn(1:3, 1));

% test that output is a string
assert(ischar(outstr))

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function [J, grad] = lrCostFunction(theta, X, y, lambda)
%LRCOSTFUNCTION Compute cost and gradient for logistic regression with
%regularization
%   J = LRCOSTFUNCTION(theta, X, y, lambda) computes the cost of using
%   theta as the parameter for regularized logistic regression and the
%   gradient of the cost w.r.t. to the parameters.

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

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Hint: The computation of the cost function and gradients can be
%       efficiently vectorized. For example, consider the computation
%
%           sigmoid(X * theta)
%
%       Each row of the resulting matrix will contain the value of the
%       prediction for that example. You can make use of this to vectorize
%       the cost function and gradient computations.
%
% Hint: When computing the gradient of the regularized cost function,
%       there're many possible vectorized solutions, but one solution
%       looks like:
%           grad = (unregularized gradient for logistic regression)
%           temp = theta;
%           temp(1) = 0;   % because we don't add anything for j = 0
%           grad = grad + YOUR_CODE_HERE (using the temp variable)
%

h = sigmoid(X * theta);
J = sum(-y .* log(h) - (1 - y) .* log(1 - h)) / m;
% regularization for cost function
J = J + (lambda / (2 * m)) * sum(theta(2:end) .^ 2);
grad = (1 / m) * X' * (h - y);
% regularization for gradient values
grad = grad + (lambda / m) * [0; theta(2:end)];

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

grad = grad(:);

end

author: zsiciarz

repo: ml-coursera


text:
function [C] = aggreg_p(op1,n,sc);
% PURPOSE: Temporal aggregation matrix preserving input dimension
% ------------------------------------------------------------
% SYNTAX: [C] = aggreg_p(op1,n,sc);
% ------------------------------------------------------------
% OUTPUT: C: nxn temporal aggregation matrix
% ------------------------------------------------------------
% INPUT:  op1: type of temporal aggregation
%         op1=1 ---> sum (flow)
%         op1=2 ---> average (index)
%         op1=3 ---> last element (stock) ---> interpolation
%         op1=4 ---> first element (stock) ---> interpolation
%         n: number of high frequency data
%         sc: number of high frequency data points
%            for each low frequency data points (freq. conversion)
% ------------------------------------------------------------
% LIBRARY: aggreg_v
% ------------------------------------------------------------
% SEE ALSO: aggreg, temporal_agg
% ------------------------------------------------------------
% NOTE: Use aggreg_v_X for extended interpolation

% written by:
%  Enrique M. Quilis
%  Macroeconomic Research Department
%  Ministry of Economy and Competitiveness
%  <enrique.quilis@mineco.es>

% Version 1.0 [October 2010]

[op1 n sc]

% ------------------------------------------------------------
% Computing implicit numer of low-frequency data
N = fix(n/sc);

% ------------------------------------------------------------
% Generation of temporal aggregation matrix
c = aggreg_v(op1,sc);
C1 = [zeros(sc-1,sc) ; c];
C = kron(eye(N),C1);

author: Sable

repo: mcbench-benchmarks


text:
function toolboxGenDoc
% Generate documentation, must run from dir toolbox.
%
% 1) Make sure to update and run toolboxUpdateHeader.m
% 2) Update history.txt appropriately, including w current version
% 3) Update overview.html file with the version/data/link to zip:
%     edit external/m2html/templates/frame-piotr/overview.html
%
% USAGE
%  toolboxGenDoc
%
% INPUTS
%
% OUTPUTS
%
% EXAMPLE
%
% See also
%
% Piotr's Image&Video Toolbox      Version 3.20
% Copyright 2013 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

% Requires external/m2html to be in path.
cd(fileparts(mfilename('fullpath'))); cd('../');
addpath([pwd '/external/m2html']);

% delete temporary files that should not be part of release
fs={'pngreadc','pngwritec','rjpg8c','wjpg8c'};
for i=1:length(fs), delete(['images/private/' fs{i} '.*']); end

% delete old doc and run m2html
if(exist('doc/','dir')), rmdir('doc/','s'); end
dirs={'channels','classify','detector','images','filters','matlab'};
m2html('mfiles',dirs,'htmldir','doc','recursive','on','source','off',...
  'template','frame-piotr','index','menu','global','on');

% copy custom menu.html and history file
sDir='external/m2html/templates/';
copyfile([sDir 'menu-for-frame-piotr.html'],'doc/menu.html');
copyfile('external/history.txt','doc/history.txt');

% remove links to private/ in the menu.html files and remove private/ dirs
for i=1:length(dirs)
  name = ['doc/' dirs{i} '/menu.html'];
  fid=fopen(name,'r'); c=fread(fid,'*char')'; fclose(fid);
  c=regexprep(c,'
([^<]*[<]?[^<]*)private([^<]*[<]?[^<]*)
','');
  fid=fopen(name,'w'); fwrite(fid,c); fclose(fid);
  name = ['doc/' dirs{i} '/private/'];
  if(exist(name,'dir')), rmdir(name,'s'); end
end

% postprocess each html file
for d=1:length(dirs)
  fs=dir(['doc/' dirs{d} '/*.html']); fs={fs.name};
  for j=1:length(fs), postProcess(['doc/' dirs{d} '/' fs{j}]); end
end

end

function postProcess( fName )
lines=readFile(fName);
assert(strcmp(lines{end-1},'') && strcmp(lines{end},''));
% remove m2html datestamp (if present)
assert(strcmp(lines{end-2}(1:22),'
<address>Generated'));
if( strcmp(lines{end-2}(1:25),'
<address>Generated on'))
  lines{end-2}=regexprep(lines{end-2}, ...
    '
<address>Generated on .* by','
<address>Generated by');
end
% insert Google Analytics snippet to end of file
ga={ '';
  '';
  '';
  '';
  '';
  '' };
lines = [lines(1:end-3); ga; lines(end-2:end)];
% write file
writeFile( fName, lines );
end

function lines = readFile( fName )
fid = fopen( fName, 'rt' ); assert(fid~=-1);
lines=cell(10000,1); n=0;
while( 1 )
  n=n+1; lines{n}=fgetl(fid);
  if( ~ischar(lines{n}) ), break; end
end
fclose(fid); n=n-1; lines=lines(1:n);
end

function writeFile( fName, lines )
fid = fopen( fName, 'wt' );
for i=1:length(lines); fprintf( fid, '%s\n', lines{i} ); end
fclose(fid);
end

author: zouchuhang

repo: LayoutNet


text:
function[varargout]=ncload(varargin)
%NCLOAD  Load all variables from a NetCDF file and convert trajectories to cells.
%
%   NCLOAD(FILENAME) loads all variables from the NetCDF file FILENAME and
%   places them into a structure of the same name in the calling workspace.
%   The extension '.nc' in FILENAME is optional.
%
%   If FILENAME contains a full pathname, only the final portion after the
%   last '/' or '\' is used for the structure name.
%
%   Note that any hyphens '-' in the filename are replaced with underscores
%   in the structure name, as the former are not allowed in variable names.
%
%   For example, if FILENAME='/Home/data' contains variables 'num', 'lat',
%   and 'lon', then the result of calling NCLOAD('data') will be a
%   structure 'data' with fields 'data.num', 'data.lat', and 'data.lon'.
%
%   Then USE can be used to map these variables into the main workspace.
%
%   NCLOAD(FILENAME,VAR1,VAR2,...,VARN) will only load the variables with
%   the names VAR1, VAR2,...VARN.
%
%   NCLOAD will convert any variables that are numeric but not doubles to
%   doubles.  This is because experience has shown using other data types
%   in code that is expecting doubles can lead to errors in Matlab.
%   __________________________________________________________________
%
%   Convert trajectory data
%
%   In addition, if FILENAME has the global attribute
%
%         featureType    = 'trajectory'
%
%   then any variables in FILENAME are have 'Dimensions: obs' will be
%   interpreted as concatentated trajectory data, following the NetCDF
%   Climate and Forecast (CF) Metadata Conventions Appendix H.4.3
%
%         http://cfconventions.org/Data/cf-conventions/cf-conventions-1.8/cf-conventions.html
%
%   NCLOAD adopts the convention that a variable named 'ids', having
%   'Dimensions: obs', provides the id for all other variables.
%
%   All variables of dimension obs will be converted to cell arrays by
%   first calling COLBREAKS using ids as the first argument, then COL2CELL.
%
%   Alternatively, all variables with 'Dimensions: obs'  will be
%   interpreted as being in the NaN-separated column format used by
%   COL2CELL, which will be called to put then in cell array format. This
%   format is supported to provide for reverse compatibility; it is
%   recommended to use the CF trajectory convention instead.
%   __________________________________________________________________
%
%   The following formats also work:
%
%          ncload filename
%          ncload filename var1 var2 ... varN
%
%   In this format, the input strings have to be the actual names of the
%   file and variables, as opposed to variables containing those names.
%
%   This is basically a way to conveniently work with small NetCDF files as
%   if they were mat-files.  Typically, NetCDF files load much faster.
%
%   Usage: ncload(filename);
%          ncload filename
%          ncload(filename,'num','lat','lon');
%          ncload filename num lat lon
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2019--2020 J.M. Lilly --- type 'help jlab_license' for details

filename=varargin{1};
structname=filename;
ii=findstr(structname,'.');
if isempty(ii)
    filename=[filename '.nc'];
else
    structname=structname(1:ii-1);
end
index=findstr(structname,'/');
if ~isempty(index)
    structname=structname(index(end)+1:end);
end
index=findstr(structname,'\');
if ~isempty(index)
    structname=structname(index(end)+1:end);
end
structname(real(structname)==real('-'))='_';  %Replace minus signs with underscores

vars=varargin(2:end);
infostruct=ncinfo(filename);

%infostruct.Dimensions

for i=1:length(infostruct.Variables)
    allvars{i}=infostruct.Variables(i).Name;
    alldimensions{i}=infostruct.Variables(i).Dimensions.Name;
end

%alldimensions

dimensions=[];
if isempty(vars)
    vars=allvars;
    dimensions=alldimensions;
else
    for i=1:length(vars)
        for j=1:length(allvars)
            if strcmp(vars{i},allvars{j})
                 dimensions{i}=alldimensions{j};
            end
        end
    end
end

try
    attvalue=ncreadatt(filename,'/','featureType');
catch
    attvalue='nontrajectory';
end

convert=false;
if strcmp(attvalue,'trajectory')
    convert=true;
    isobs=false(length(dimensions),1);
    isids=false(length(dimensions),1);
    for i=1:length(isobs)
        isobs(i)=strcmp(dimensions{i},'obs');
        isids(i)=strcmp(vars{i},'ids');
    end
    if anyany(isobs)
        if ~anyany(isids)
            convert=false;
            disp('NCLOAD finding trajectory data but no "ids" variable, so not converting.')
        end
    end

end

if ~isempty(dimensions)
    for i=1:length(vars)
        if ~isempty(dimensions{i})
            %infostruct.Variables(i).Dimensions.Name
            if strcmp(dimensions{i},'column')
                %converting my older 'column' format
                str=[structname '.' vars{i} '=doubifnum_internal(col2cell_internal(ncread(''' filename ''',''' vars{i} ''')));'  ];
            else
                str=[structname '.' vars{i} '=doubifnum_internal(ncread(''' filename ''',''' vars{i} '''));'  ];
            end
            evalin('caller',str)
        end
    end
end

if convert
    strlist=[structname '.ids'];%make sure ids comes firsts
    for i=1:length(vars)
        if strcmp(dimensions{i},'obs')&&~strcmp(vars{i},'ids')
            strlist=[strlist ',' structname '.' vars{i} ];
        end
    end
    evalin('caller',['[' strlist ']=colbreaks_internal(' strlist ');'])
    evalin('caller',['[' strlist ']=col2cell_internal(' strlist ');'])
end

function[x]=doubifnum_internal(x)
%DOUBIFNUM  Converts a variable to double precision if it is numeric.
%
%   X=DOUBIFNUM
%
%   Usage: x=doubifnum(x);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2020 J.M. Lilly --- type 'help jlab_license' for details

if isnumeric(x)
    x=double(x);
end

function[varargout]=colbreaks_internal(varargin)
%COLBREAKS  Insert NANs into discontinuties in a vector.
%
%   O1=COLBREAKS(I1), where I1 is a column vector, inserts NANs into I1 at
%   its discontinuities (and puts one at the end), for example:
%
%	  COLBREAKS([2;2;2;3;3])=[2;2;2;NAN;3;3;NAN];
%
%   [O1,O2,O3,...]=COLBREAKS(I1,I2,...), where I1,I2,... are all column
%   vectors of the same length, uses I1 as a reference for the other
%   vectors and inserts NANs into all the vectors at locations where I1 is
%   discontinous.
%
%   For instance, using station number for I1 will sort a column of CTD
%   data into a NAN-padded matrix.
%
%   Note that for complex-valued input arrays, a complex-valued NAN,
%   NAN+SQRT(-1)*NAN, is used to indicate the break locations.
%
%   MAT2=COLBREAKS(MAT1), where MAT1 and MAT2 are matrices of the same
%   size, also works.  In this case the first column of MAT1 is used as the
%   reference vector.
%
%   COLBREAKS, COL2MAT, and MAT2COL together form a system for moving data
%   with segments of nonuniform length rapidly back and forth between a
%   column format and a padded-matrix format. CTD or float data, for
%   instance, can be stored in the (usually much smaller) column format and
%   converted into the matrix format upon loading.
%
%   COLBREAKS(C1,C2,...); with no output arguments overwrites the original
%   input variables.
%
%   See also COL2MAT, MAT2COL, COL2CELL, CELL2COL, ORBITBREAKS.
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2000--2013 J.M. Lilly --- type 'help jlab_license' for details

bmat=0;
nargs=nargin;
ref=varargin{1};
if nargin==1 && size(ref,2)>1
    temp=ref;
    ref=temp(:,1);
    c=cell(size(temp,2),1);
    for i=1:size(temp,2)-1
        c{i}=temp(:,i+1);
    end
    nargs=size(temp,2);
    clear temp
    bmat=1;
else
    c=varargin(2:end);
end

ii=0;
index=find(diff(ref)~=0)+1;
bool=false(size(ref));
bool(index)=ones(size(index));
index=(1:length(ref))'+cumsum(bool);

while ii<nargs - 1
    ii=ii+1;
    col=c{ii};

    if isreal(col)
        colout=nan*ones(max(index),1);
        colout(index)=col;
        colout=[colout;nan];
    else
        colout=(nan+sqrt(-1)*nan)*ones(max(index),1);
        colout(index)=col;
        colout=[colout;(nan+sqrt(-1)*nan)];
    end

    varargout{ii+1}=colout;
end

if isreal(ref)
    refout=nan*ones(max(index),1);
    refout(index)=ref;
    refout=[refout;nan];
else
    refout=(nan+sqrt(-1)*nan)*ones(max(index),1);
    refout(index)=ref;
    refout=[refout;(nan+sqrt(-1)*nan)];
end

varargout{1}=refout;
if bmat
    for i=1:nargs-1
        varargout{1}(:,i+1)=varargout{i+1};
    end
    varargout=varargout(1);
end

function[varargout]=col2cell_internal(varargin)
%COL2CELL  Converts 'column-appended' data into cell arrays of numeric arrays.
%
%   X=COL2CELL(COL) converts the array COL, a 2D array having blocks of
%   data separated by rows of all NaNs, into a cell array, with each cell
%   containing one of the data blocks, excluding the trailing NaNs.
%
%   Since NANs mark the end of each block, missing or bad data in COL
%   should instead be indicated by INFs.
%
%   [X1,X2,...,XN]=COL2CELL(C1,C2,...,CN) also works, where the input
%   fields C1,...,CN are all the same size.  In this case the locations of
%   NANs in C1 are used to break all the other CN into cells.
%
%   COL2CELL(C1,C2,...,CN); with no output arguments overwrites the
%   original variables.
%
%   COL2CELL is inverted by CELL2COL.
%   __________________________________________________________________
%
%   See also CELL2COL, COL2MAT, MAT2COL, COLBREAKS.
%
%   'col2cell --t' runs a test.
%
%   Usage: x=col2cell(col);
%          [x1,x2,...,xN]=col2cell(c1,c2,...,cN);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2008--2018 J.M. Lilly --- type 'help jlab_license' for details

i=0;
if ~isempty(varargin{1})
    while i<nargin
        i=i+1;
        col=varargin{i};
        if i==1
            %Note minor changes to make this work on matrices rather than
            %just columns.
            ia=[1;find(isnan(col(1:end-1,1)))+1];
            ib=find(isnan(col(:,1)))-1;
            %[L,ia,ib]=blocklen(col);
            ray=cell(length(ia),1);
        end
        for j=1:length(ia)
            if ~isempty(col)
                ray{j}=col(ia(j):ib(j),:);
            else
                ray{j}=[];
            end
        end
        varargout{i}=ray;
    end
else
    for i=1:nargout
        varargout{i}=varargin{1};
    end
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function prob = rowlin2nl(prob,sparsity,warn)
%ROWLIN2NL Convert Row Based Linear Constraints to Row Based Nonlinear Constraints

%   Copyright (C) 2012 Jonathan Currie (IPL)

if(isempty(prob.A))
    return; %nothing to do
end

if(~exist('warn','var') || isempty(warn))
    warn = 1;
end
if(~exist('sparsity','var') || isempty(sparsity))
    sparsity = 1;
end

%If we don't have any nonlinear constraints, just build from linear ones
if(isempty(prob.nlcon))
    %Build nonlinear constraints + jacobian + jac structure (assumes jac + str not provided)
    prob.nlcon = @(x) prob.A*x;
    if(sparsity)
        prob.nljac = @(x) sparse(prob.A);
        prob.nljacstr = @() sparse(double(prob.A ~= 0));
    else
        prob.nljac = @(x) prob.A;
        prob.nljacstr = @() double(prob.A ~= 0);
    end
    %Assign constraint bounds
    prob.cl = prob.rl;
    prob.cu = prob.ru;
    %Add sizes
    eq = prob.cl == prob.cu; neq = ~eq;
    ile = isfinite(prob.cu) & neq;
    ige = isfinite(prob.cl) & neq;
    prob.sizes.nnlineq = sum(ile + ige);
    prob.sizes.nnleq = sum(eq);
    %Indicate all linear
    prob.lincon = 1;

%Otherwise Augment linear to nonlinear constraints
else
    %Augment nonlinear constraint function
    prob.nlcon = @(x) [prob.nlcon(x);
                       prob.A*x];
    %Augment nonlinear Jacobian if it exists
    if(~isempty(prob.nljac))
        if(sparsity)
            A = sparse(prob.A);
            prob.nljac = @(x) [prob.nljac(x);
                               A];
        else
            prob.nljac = @(x) [prob.nljac(x);
                               prob.A];
        end
    end
    %Augment nonlinear jacobian structure if it exists
    if(~isempty(prob.nljacstr))
        if(sparsity)
            A = sparse(prob.A ~= 0);
            prob.nljacstr = @() [prob.nljacstr(); A];
        else
            prob.nljacstr = @() [prob.nljacstr(); double(prob.A ~= 0)];
        end
    end

    %Augment rhs + type
    prob.cl = [prob.cl;
               prob.rl];
    prob.cu = [prob.cu;
               prob.ru];
    %Add sizes
    eq = prob.rl == prob.ru; neq = ~eq;
    ile = isfinite(prob.ru) & neq;
    ige = isfinite(prob.rl) & neq;
    prob.sizes.nnlineq = prob.sizes.nnlineq + sum(ile + ige);
    prob.sizes.nnleq = prob.sizes.nnleq + sum(eq);
    %Indicate not all linear
    prob.lincon = 0;
end

if(warn > 1)
    optiwarn('opti:rowlin_nl','RowLin2NL - Linear Constraints have been converted to Nonlinear Constraints to suit NLP Solver interface');
end

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function geometry_test039 ( )

%*****************************************************************************80
%
%% TEST039 tests LINES_EXP_INT_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST039\n' );
  fprintf ( 1, '  LINES_EXP_INT_2D finds intersections of\n' );
  fprintf ( 1, '    two explicit lines in 2D.\n' );
  fprintf ( 1, '\n' );

  for itest = 1 : 3
%
%  x + 2y - 4 = 0
%  x - y - 1 = 0
%
    if ( itest == 1 )

      p1(1:2,1) = [ 0.0;  2.0 ];
      p2(1:2,1) = [ 4.0;  0.0 ];
      q1(1:2,1) = [ 0.0; -1.0 ];
      q2(1:2,1) = [ 1.0;  0.0 ];
%
%  x + 2y - 4 = 0
%  2x + 4y - 1 = 0
%
    elseif ( itest == 2 )

      p1(1:2,1) = [ 0.00; 2.00 ];
      p2(1:2,1) = [ 4.00; 0.00 ];
      q1(1:2,1) = [ 0.00; 0.25 ];
      q2(1:2,1) = [ 0.50; 0.00 ];
%
%  x + 2y - 4 = 0
%  -3x - 6y +12 = 0
%
    elseif ( itest == 3 )

      p1(1:2,1) = [ 0.0; 2.0 ];
      p2(1:2,1) = [ 4.0; 0.0 ];
      q1(1:2,1) = [ 0.0; 2.0 ];
      q2(1:2,1) = [ 4.0; 0.0 ];

    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  P1  %8f  %8f\n', p1(1:2,1) );
    fprintf ( 1, '  P2  %8f  %8f\n', p2(1:2,1) );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Q1  %8f  %8f\n', q1(1:2,1) );
    fprintf ( 1, '  Q2  %8f  %8f\n', q2(1:2,1) );

    [ ival, p ] = lines_exp_int_2d ( p1, p2, q1, q2 );

    if ( ival == 1 )
      fprintf ( 1, '  Intersection at %8f  %8f\n', p(1:2,1) );
    elseif ( ival == 0 )
      fprintf ( 1, '  Lines are parallel, no intersection.\n' );
    elseif ( ival == 2 )
      fprintf ( 1, '  Lines are coincident.\n' );
    else
      fprintf ( 1, '  Unknown return value of IVAL = %d\n', ival );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function scn_export_papersetup(minsize)
% :Usage:
% ::
%
%    scn_export_papersetup([opt: min size in pixels, default = 400])
%
% set paper size for current figure so that print to png or tiff looks as
% it should (as it does on-screen)
%
% ..
%    tor wager, aug. 06
% ..

if nargin < 1, minsize = 400; end

% make sure that min size of fig. is 400 pixels
sz = get(gcf,'Position');
sz = sz(3:4);   % width and height
szratio = max(sz) ./ min(sz);   % max to min size
wh = find(sz == min(sz));
sz(wh) = minsize;
wh = find(sz == max(sz));
sz(wh) = minsize * szratio;

% set paper size using current screen
set(gcf,'PaperUnits','points','PaperPosition',[0 0 sz],'PaperType','usletter');
return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function newBb = convertBb_X1Y1X2Y2_to_X1Y1WH( bb )
%convertBb_X1Y1X2Y2_to_X1Y1WH converts the bounding boxes from [X1,Y1,X2,Y2] format to [X1,Y1,W,H] format

numBb = size(bb, 1);

newBb = zeros(numBb, 4);

newBb(:, 1) = bb(:, 1);
newBb(:, 2) = bb(:, 2);
newBb(:, 3) = bb(:, 3) - bb(:, 1) + 1;
newBb(:, 4) = bb(:, 4) - bb(:, 2) + 1;

end

author: aosokin


text:
function [ vX ] = HybridOrthogonalProjectionOntoConvexSets( cProjFun, vY, numIterations, stopThr, compMethod )
% ----------------------------------------------------------------------------------------------- %
% [ vX ] = HybridOrthogonalProjectionOntoConvexSets( cProjFun, vY, numIterations, stopThr, compMethod )
%   Solves \arg \min_{x} 0.5 || x - y ||, s.t. x \in \bigcap {C}_{i} using
%   Hybrid Orthogonal method.
% Input:
%   - cProjFun      -   Array of Projection Functions.
%                       Cell array of anonymouse functions which each is a
%                       projection into a sub space.
%                       Structure: Cell Array.
%                       Type: NA.
%                       Range: NA.
%   - vY            -   Input Vector.
%                       Input vector to be projected.
%                       Structure: Vector (m x 1).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
%   - numIterations -   Number of Iterations.
%                       Sets the number of iterations of the algorithm.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: {1, 2, ...}.
%   - stopThr       -   Stopping Threshold.
%                       Sets the threshold between consecutive iterations
%                       for stopping the algorithm.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: (0, inf).
%   - compMethod    -   Hybrid Method.
%                       Sets the hybrid method to calculate.
%                       Structure: Scalar.
%                       Type: 'Single' / 'Double'.
%                       Range: {1, 2}.
% Output:
%   - vX            -   Solution Vector.
%                       The solution to the optimization problem.
%                       Structure: Vector (m x 1).
%                       Type: 'Single' / 'Double'.
%                       Range: (-inf, inf).
% References
%   1.  Quadratic Optimization of Fixed Points of Non Expensive Mappings in Hilbert Space (https://doi.org/10.1080/01630569808816822).
% Remarks:
%   1.  B
% TODO:
%   1.  C
% Release Notes:
%   -   1.0.000     19/03/2020  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

COMP_METHOD_A = 1;
COMP_METHOD_B = 2;

switch(compMethod)
    case(COMP_METHOD_A)
        vX = HybridOrthogonalProjectionOntoConvexSetsA(cProjFun, vY, numIterations, stopThr);
    case(COMP_METHOD_B)
        vX = HybridOrthogonalProjectionOntoConvexSetsB(cProjFun, vY, numIterations, stopThr);
end

end

function [ vX ] = HybridOrthogonalProjectionOntoConvexSetsA( cProjFun, vY, numIterations, stopThr )

numSets     = size(cProjFun, 1);
numElements = size(vY, 1);

vX = vY;
vU = vX;

kk = 0;

for ii = 1:numIterations

    vU(:) = vX;

    % The loop is as vX changes slowly or not at all with each projection.
    % Hence in order to prevent pre mature stopping one should apply all
    % projections and then generate new vX.
    for jj = 1:numSets
        kk = kk + 1;
        % See Quadratic Optimization of Fixed Points of Non Expensive
        % Mappings in Hilbert Space (Page 12, Equation 20, The note after Equation 21)
        % paramLambdaN = 1 / ((kk + 1) ^ 0.99);
        % See Quadratic Optimization of Fixed Points of Non Expensive
        % Mappings in Hilbert Space (Page 14, Equation 30)
        % Pay attention that Equation 30 allows use of 1 / n.
        paramLambdaN = 1 / (kk + 1);
        % This still require the Fixed Point of the intersection is the
        % same for any cyclic variation of the Sets (Which happens for nay
        % Non Expensive Projection).
        vX(:) = (paramLambdaN * vY) + (1 - paramLambdaN) * cProjFun{jj}(vX);
    end

    stopCond = sum(abs(vU - vX)) < stopThr;

    if(stopCond)
        break;
    end
end

end

function [ vX ] = HybridOrthogonalProjectionOntoConvexSetsB( cProjFun, vY, numIterations, stopThr )

numSets     = size(cProjFun, 1);
numElements = size(vY, 1);

vX = vY;
vU = vX;
vT = vX;

vW = rand(numSets, 1);
vW = vW / sum(vW);

for ii = 1:numIterations

    vU(:) = vX;

    % See Quadratic Optimization of Fixed Points of Non Expensive Mappings
    % in Hilbert Space (Page 18, Equation 44)
    vT(:) = 0;
    for jj = 1:numSets
        vT(:) = vT + (vW(jj) * cProjFun{jj}(vX));
    end
    paramLambdaN = 1 / (ii + 1);
    vX(:) = (paramLambdaN * vY) + (1 - paramLambdaN) * vT;

    stopCond = sum(abs(vU - vX)) < stopThr;

    if(stopCond)
        break;
    end
end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
clear
close all

plotx2east

load martensite_single_grain

[grains,ebsd('indexed').grainId] =...
    calcGrains(ebsd('indexed'),'angle',3*degree);

cond_grains = ismember(grains.id,unique(ebsd(ids_of_interest).grainId));

%%
job = parentGrainReconstructor(ebsd,grains) % define job
job.useBoundaryOrientations = true;
clearvars grains ebsd

job.p2c = orientation.KurdjumovSachs(job.p2c.CS,job.p2c.SS);
job.calcParent2Child('quantile',0.6,'c2c'); % get representative OR

%%
job.calcHyperGraph3('threshold',5*degree,'c2c','mergeSimilar','mergethreshold',8*degree,'keepGraph')

numIters = [0 1 3 10 50];
grainId = 4691;

%%

for k = 1:length(numIters);
    if k == 1
    job.clusterHyperGraph3('numIter',numIters(k),...
        'inflationPower',1,'merged','mergethreshold',8*degree,'keepGraph')
        graphplotfunc(job,cond_grains,grainId,k,numIters)
    else
    job.clusterHyperGraph3('numIter',numIters(k) - numIters(k-1),...
        'inflationPower',1,'merged','mergethreshold',8*degree,'keepGraph')
        graphplotfunc(job,cond_grains,grainId,k,numIters)
    end
end

%%
p2c = orientation.NishiyamaWassermann(job.p2c.CS,job.p2c.SS);
p2c = p2c.project2FundamentalRegion(job.p2c);
job.p2c = p2c;

job.calcParentFromVote;
nextAxis
plot(job.grains(cond_grains & job.grains.phase == 2),...
    job.grains(cond_grains & job.grains.phase == 2).meanOrientation)
lims = job.ebsd.extent;
xlim([lims(1) lims(2)])
ylim([lims(3) lims(4)])

function graphplotfunc(job,cond_grains,grainId,k,numIters)

p2c = orientation.NishiyamaWassermann(job.p2c.CS,job.p2c.SS);
p2c = p2c.project2FundamentalRegion(job.p2c);

numV = length(variants(p2c,'parent'));

[i,j] = find(tril(job.graph));
i1 = ceil(i/numV);
i2 = i-(i1-1)*numV;
j1 = ceil(j/numV);
j2 = j-(j1-1)*numV;
prob = nonzeros(tril(job.graph));

cond = i1 ~= j1;

i1 = i1(cond);
j1 = j1(cond);

cond2 = (i1 == grainId | j1 == grainId);
i1 = i1(cond2);
j1 = j1(cond2);

i2 = i2(cond);
j2 = j2(cond);

i2 = i2(cond2);
j2 = j2(cond2);

prob = prob(cond);
prob = prob(cond2);

pOri = variants(p2c, job.grains(j1).meanOrientation, j2);
ipfKey = ipfHSVKey(job.csParent);
colors = ipfKey.orientation2color(pOri);

%%
%Equivalency:
z = unique([i1;j1]);
ord = zeros(max(z),1);
ord(z) = 1:length(z);

[x,y] = centroid(job.grains);
G = graph(ord(i1),ord(j1),prob);

if k == 1
    newMtexFigure('layout',[2,3])
else
    nextAxis
end

plot(job.grains(cond_grains).boundary)
hold on
plot(job.grains(unique([i1;j1])).boundary,'linecolor','red')
hold on
plot(G,'XData',x(unique([i1;j1])),'YData',y(unique([i1;j1])),...
    'EdgeColor',colors,'NodeColor','none','NodeLabel',{})
hold on
lims = job.ebsd.extent;
text(lims(2)-40,lims(3)+15, {[num2str(numIters(k)) ' iterations'],...
    [num2str(length(z)) ' nodes'],...
    [num2str(length(i1)) ' edges']}, 'HorizontalAlignment','left',...
    'FontSize',20, 'FontWeight','bold','BackgroundColor', 'white')

xlim([lims(1) lims(2)])
ylim([lims(3) lims(4)])

end

author: mtex-toolbox

repo: mtex

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