text:
function jed = ymdf_to_jed_alexandrian ( y, m, d, f )

%*****************************************************************************80
%
%% YMDF_TO_JED_ALEXANDRIAN converts an Alexandrian YMDF date to a JED.
%
%  Discussion:
%
%    This code needs to be adjusted to fit the Alexandrian model.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, M, D, real F, the YMDF date.
%
%    Output, real JED, the corresponding Julian Ephemeris Date.
%

%
%  Convert the calendar date to a computational date.
%
  y_prime = y + 4690 - floor ( ( 13 - m ) / 13 );
  m_prime = mod ( m + 12, 13 );
  d_prime = d - 1;
%
%  Convert the computational date to a JED.
%
  jed = floor ( ( 1461 * y_prime ) / 4 ) + 30 * m_prime + d_prime - 124 - 0.5;
  jed = jed + f;

  return
end


text:
%  Figure 4.6      Feedback Control of Dynamic Systems, 6e
%                   Franklin, Powell, Emami
%
numG=2;
den=[1 1 2];
% add integrator
denG=conv(den,[1 0]);
t=0:.03:10;
y=step(numG,denG,t);
axis('square');
plot(t,t,t,y);
title('Fig. 4.6 Relationship between ramp response and Kv');
xlabel('Time (sec)');
ylabel('r, y');
axis('normal')
gtext('r')
gtext('y')
gtext('e_{ss}=1/K_v')
nicegrid;

author: Sable

repo: mcbench-benchmarks


text:
function [density, pressure, temperature] = getAtmoDensityAtAltitude(bodyInfo, altitude, lat, ut, long)
%getAtmoDensityAtAltitude Summary of this function goes here
%   Detailed explanation goes here
%   See here for math: https://forum.kerbalspaceprogram.com/index.php?/topic/142686-modeling-atmospheres-in-ksp/

    if(altitude <= bodyInfo.atmohgt && altitude >= 0)
%         pressure = getPressureAtAltitude(bodyInfo, altitude);
        pressure = bodyInfo.getBodyAtmoPressure(altitude);

        if(pressure > 0)
            temperature = bodyInfo.getBodyAtmoTemperature(ut, lat, long, altitude);

            density = getDensityFromIdealGasLaw(pressure, temperature, bodyInfo.atmomolarmass);
            if(density < 0)
                density = 0;
                pressure = 0;
                temperature = 0;
            end
        else
            density = 0;
            pressure = 0;
            temperature = 0;
        end
    elseif(altitude <= 0)
        density = 0;
        pressure = 0;
        temperature = 0;
    else
        density = 0;
        pressure = 0;
        temperature = 0;
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function determ = symm_random_determinant ( n, d, key )

%*****************************************************************************80
%
%% SYMM_RANDOM_DETERMINANT returns the determinant of the SYMM_RANDOM matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 March 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real D(N), the desired eigenvalues for the matrix.
%
%    Input, integer KEY, a positive integer that selects the data.
%
%    Output, real DETERM, the determinant.
%
  determ = prod ( d(1:n) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function ppm = ParforProgressStarter2(s, n, percentage, do_debug)
% Starter function for parfor progress-monitor, that will automatically
% choose between a text or a GUI progress monitor. You can use this
% progress-monitor in both for and parfor loops.
%
% ParforProgressStarter2 is the successor of ParforProgressStarter - it is
% roughly 12 times more performant during GUI mode!
%
% you should use it like this:
%
%   N = 1000;
%   try % Initialization
%       ppm = ParforProgressStarter2('test', N, 0.1);
%   catch me % make sure "ParforProgressStarter2" didn't get moved to a different directory
%       if strcmp(me.message, 'Undefined function or method ''ParforProgressStarter2'' for input arguments of type ''char''.')
%           error('ParforProgressStarter2 not in path.');
%       else
%           % this should NEVER EVER happen.
%           msg{1} = 'Unknown error while initializing "ParforProgressStarter2":';
%           msg{2} = me.message;
%           print_error_red(msg);
%           % backup solution so that we can still continue.
%           ppm.increment = nan(1, nbr_files);
%       end
%   end
%
%   parfor i = 1 : N
%       your_crazy_function();
%       ppm.increment(i);
%   end
%
%   try % use try / catch here, since delete(struct) will raise an error.
%       delete(ppm);
%   catch me %#ok<NASGU>
%   end
%
%
% Copyright (c) 2010-2012, Andreas Kotowicz
%
    %%

    if nargin < 2
        disp('usage: ppm = ParforProgressStarter2( text, number_runs, update_percentage)');
        ppm = [];
        return;
    end

    if nargin < 3
        percentage = 0.1;
    end

    if nargin < 4
        do_debug = 0;
    end

    %% determine whether java and awt are available
    java_enabled = 1;
    if usejava('jvm') == 0
        java_enabled = 0;
    end

    awt_available = 1;
    if usejava('awt') == 0
        awt_available = 0;
    end

    % is matlab pool active?
    pool_slaves = pool_size();

    %% check for different usage scenarios:

    % 1 = ParforProgress2(GUI)
    % 2 = ParforProgress2(CONSOLE)
    % 3 = fallback ParforProgressConsole2

    % by default we use the GUI
    version_to_use = 1;

    if java_enabled == 0 % no java -> use console only
        version_to_use = 3;
    else
        if awt_available == 0 % java but no AWT -> java console
            version_to_use = 2;
        end
    end

    % check for old matlab version
    % - everything before 2008b (7.7) can cause trouble (because of saveobj /
    % loadobj)
    if get_matlab_version() < 7.07
        version_to_use = 3;
    end

    %% add directory to javapath and path
    a = which(mfilename);
    dir_to_add = fileparts(a);

    if pool_slaves > 0
        if java_enabled == 1
            pctRunOnAll(['javaaddpath({''' dir_to_add '''})']);
        end
        pctRunOnAll(['addpath(''' dir_to_add ''')']);
    else
        if java_enabled == 1
            javaaddpath({dir_to_add});
        end
        addpath(dir_to_add);
    end

    %%
    switch version_to_use
        case 1
            use_gui = 1;
            ppm = ParforProgress2(s, n, percentage, do_debug, use_gui);
        case 2
            use_gui = 0;
            ppm = ParforProgress2(s, n, percentage, do_debug, use_gui);
        case 3
            % no java, no awt, or old matlab version
            disp('Progress will update in arbitrary order.');
            ppm = ParforProgressConsole2(s, n, percentage);
        otherwise
            disp('not defined');
            ppm = [];
    end

end
%% EOF

author: thomaskuestner


text:
%% An example run of VBRPCA
clear
clc

%% Create a low-rank + sparse matrix
% Dimensions
m = 500;
n = 500;
r = 25; % rank of the low-rank component

A = randn(m,r);
B = randn(r,n);
X_true = A*B; % low rank component

SparseRatio = 0.05;
E_true = zeros(size(X_true));
Eomega = randsample(m*n, round(m*n*SparseRatio));
E_true(Eomega) = -10+20*rand(length(Eomega),1); % sparse component

% Observation
Y = X_true + E_true;

% Add noise?
sigma = 0; % no noise
% sigma = 1e-3; % some noise
% sigma = sqrt(1e-3); % more noise

Y = Y + sigma*randn(size(Y));

SNR = 10*log10( sum(sum((X_true+E_true).^2)) / sum(sum((Y-X_true-E_true).^2)));

%% Run VRBPCA
% all options are *optional*, everything will be set automatically
% you can modify these options to get better performance
options.verbose = 1;
options.initial_rank = 'auto'; % This sets to the maximum possible rank
% options.initial_rank = 300; % or we can use a value.
options.X_true = X_true;
options.E_true = E_true;
options.inf_flag = 2; % inference flag for the sparse component
% 1 for standard VB, 2 for MacKay. MacKay generally converges faster.

options.MAXITER = 200;
%Estimate noise variance? (beta is inverse noise variance)
options.UPDATE_BETA = 1;
% If the noise inv. variance is not to be estimated, set
% options.UPDATE_BETA = 0; % and set beta using
% options.beta = 1e3;

% Select the optimization mode:
% 'VB': fully Bayesian inference (default)
% 'VB_app': fully Bayesian with covariance approximation
% 'MAP': maximum a posteriori (covariance is set to 0)
options.mode = 'VB';
% For large scale problems, set to 'VB_app'.
% options.mode = 'VB_app';

fprintf('Running VBRPCA...\n')
tic
[X_hat, A_hat, B_hat, E_hat] = VBRPCA(Y,options);
t_total = toc;

% Results
Xerr = norm( X_hat - X_true, 'fro' ) / norm( X_true, 'fro');
Eerr = norm( E_hat - E_true, 'fro' ) / norm( E_true, 'fro');

r_hat = rank(X_hat);
% the model does not allow for "exact" sparsity, but only "soft" sparsity.
% That is, most coefficients are almost zero (very small values), but not
% exactly zero. If the sparse coefficients are needed, one can use
% thresholding such as
nonzero_E = length(find(abs(E_hat)>1e-7));

% Show results
fprintf('\n\n-------------------------------------------------------------------\n')
fprintf('Dimensions (%d, %d), true rank = %d, nonzero E elements = %d, SNR = %g\n',m, n, r, length(Eomega), SNR);
fprintf('low rank comp. error = %g, sparse comp. error %g, est. rank = %d,\nnonzero E elements = %d, running time = %g\n', Xerr, Eerr, r_hat, nonzero_E, t_total);
fprintf('-------------------------------------------------------------------\n')

author: andrewssobral

repo: lrslibrary


text:
function solution = MP_GWO(UAV, SearchAgents, Max_iter)
%MP_GWO 多种群灰狼优化算法
%Multi Population Gray Wolf Optimization

% 超参数
g = 50;       % 动态更新加权系数

% 算法初始化
[WolfPops, ~] = PopsInit(UAV, SearchAgents, false);   % 随机生成 初始狼群
ClassPops = PopsCluster(WolfPops, UAV);                    % 进行初始聚类（用来获得 k 值）
dim = WolfPops.PosDim;                                                 % 状态变量维度
cSearchAgents = ClassPops.SearchAgents;                    % 搜索智能体个数（子种群数量）
SearchAgents = cSearchAgents * ClassPops.k;              % 对所有智能体个数进行修正（k的整数倍）
WolfPops.Pos = WolfPops.Pos(1:SearchAgents, :);       % 对种群进行修正

% 报错
if cSearchAgents < 4
    error('搜索智能体个数过少')
end

% 初始化解
Alpha_pos = zeros(ClassPops.k, dim);          % α解
Alpha_score = 1 ./ zeros(ClassPops.k, 1);     % α解适应度

Beta_pos = zeros(ClassPops.k, dim);            % β解
Beta_score = 1 ./ zeros(ClassPops.k, 1);       % β解适应度

Delta_pos = zeros(ClassPops.k, dim);           % δ解
Delta_score = 1 ./ zeros(ClassPops.k, 1);      % δ解适应度

Fitness_list = zeros(ClassPops.k, Max_iter);  % 适应度曲线

Pops.PosDim = WolfPops.PosDim;  % 子种群
Pops.lb = WolfPops.lb;
Pops.ub = WolfPops.ub;

% 迭代求解
tic
fprintf('>>MP-GWO 优化中    00.00%%')
for iter = 1 : Max_iter

    % ①  更新参数a
    a = 2 - iter * 2 / Max_iter;               % 线性递减
    %a = 2 * cos((iter/Max_iter)*pi/2);   % 非线性递减

    % ②  聚类
    if iter > 1 %初始化时聚完一次了，节省计算量
    ClassPops = PopsCluster(WolfPops, UAV);
    end
    for k = 1 : ClassPops.k
        Positions = ClassPops.Pos{k};

        % ③  寻找 α、β、δ 狼
        for i = 1 : cSearchAgents
            % 读取目标函数
            fitness = ClassPops.Fitness(k, i);

            % 更新 Alpha、Beta 和 Delta 解
            if fitness <= Alpha_score(k)  % 适应能力最强（因为性能指标越小越好，因此为小于号）
                Alpha_score(k) = fitness;
                Alpha_pos(k, :) = Positions(i, :);
            end
            if fitness > Alpha_score(k) && fitness <= Beta_score(k)
                Beta_score(k) = fitness;
                Beta_pos(k, :) = Positions(i, :);
            end
            if fitness > Alpha_score(k) && fitness > Beta_score(k) && fitness <= Delta_score(k)
                Delta_score(k) = fitness;
                Delta_pos(k, :) = Positions(i, :);
            end
        end

        % ④  更新位置（朝着前三只狼位置前进）
        for i = 1 : cSearchAgents
            for j = 1 : dim

                r1 = rand();
                r2 = rand();
                A1 = 2*a*r1 - a;
                C1 = 2*r2;
                D_alpha = abs(C1*Alpha_pos(k, j) - Positions(i, j));
                X1 = Alpha_pos(k, j) - A1*D_alpha;

                r1 = rand();
                r2 = rand();
                A2 = 2*a*r1 - a;
                C2 = 2*r2;
                D_beta = abs(C2*Beta_pos(k, j) - Positions(i, j));
                X2 = Beta_pos(k, j) - A2*D_beta;

                r1 = rand();
                r2 = rand();
                A3 = 2*a*r1 - a;
                C3 = 2*r2;
                D_delta = abs(C3*Delta_pos(k, j) - Positions(i, j));
                X3 = Delta_pos(k, j) - A3*D_delta;

                % 静态更新
                Positions(i, j) = (X1 + X2 + X3) / 3;
                % 动态更新
                %                 q = g * a; %阈值
                %                 if abs(Alpha_score(k)-Delta_score(k)) > q
                %                     Sum_score = Alpha_score(k) + Beta_score(k) + Delta_score(k);
                %                     Positions(i, j) = (Alpha_score(k)*X1 + Beta_score(k)*X2 + Delta_score(k)*X3) / Sum_score;
                %                 else
                %                     Positions(i, j) = (X1 + X2 + X3) / 3;
                %                 end

            end
        end

    % ⑤  调整不符合要求的状态变量
    Pops.Pos = Positions;
    ProbPoints =  ClassPops.ProbPoints{k};
    [Pops, ~] = BoundAdjust(Pops, ProbPoints, UAV);

    % ⑥  存储适应度
    Fitness_list(k, iter) = Alpha_score(k);

    % ⑦  合并种群
    WolfPops.Pos(cSearchAgents*(k-1)+1:cSearchAgents*k, :) = Pops.Pos;

    end

    if iter/Max_iter*100 < 10
        fprintf('\b\b\b\b\b%.2f%%', iter/Max_iter*100)
    else
        fprintf('\b\b\b\b\b\b%.2f%%', iter/Max_iter*100)
    end
end
fprintf('\n\n>>计算完成！\n\n')
toc

% 寻找 α β δ 位置
n = 3;
A = ClassPops.Fitness;
t = findmin(A, n);
index = cSearchAgents * (t(:, 1) - 1) + t(:, 2);

real_Alpha_no = index(1);
real_Beta_no = index(2);
real_Delta_no = index(3);
Alpha_Data = ClassPops.Data{t(1, 1)}{t(1, 2)}  ;

% 输出值
solution.method = 'MP-GWO';                                 % 算法
% solution.ClassPops = ClassPops;                               % 分类信息
solution.WolfPops = WolfPops;                                % 所有解种群信息
solution.Tracks = Pops2Tracks(WolfPops, UAV);    % 所有解航迹信息
solution.Fitness_list = mean(Fitness_list, 1);            % 所有α解的平均适应度曲线
solution.Alpha_Data = Alpha_Data;                          % 真 · α 的威胁信息
solution.Alpha_no = real_Alpha_no;                          % 真 · α 的位置
solution.Beta_no = real_Beta_no;                               % 真 · β 的位置
solution.Delta_no = real_Delta_no;                            % 真 · δ 的位置

end

%% 寻找A矩阵中最小n个数的位置
function t = findmin(A, n)
    t = sort(A(:));
    [x, y] = find(A <= t(n), n);
    t = [x, y];         % 前n个最小项在矩阵A中的位置[行,列]
    B = zeros(n, 1);
    for i = 1 : n
        B(i) = A(t(i, 1), t(i, 2));
    end
    [~, index] = sort(B);
    t = t(index, :); % 前n个从小到大排序的位置
end

%% 对种群进行聚类
function [ClassPops] = PopsCluster(WolfPops, UAV)

SearchAgents = size(WolfPops.Pos, 1);  % 智能体个数
Dim = WolfPops.PosDim;                        % 智能体维度
Tracks = Pops2Tracks(WolfPops, UAV); % 智能体转换成航迹信息

% 计算适应度
o_Fitness = zeros(SearchAgents, 1); % 60*1
o_subF = []; % 5*60
o_ProbPoints = cell(SearchAgents, 1);  % 60*1
o_Data = cell(SearchAgents, 1); % 60*1

% parfor 并行计算适应度
for i = 1:SearchAgents
    [fitness, subF, Data] = ObjFun(Tracks{i}, UAV);
    o_ProbPoints{i} = Data.ProbPoint;  %cell-cell
    o_Data(i) = {Data}; %cell-struct
    o_Fitness(i) = fitness; %vector-var
    o_subF = [o_subF, subF];
end

% 分类
k = size(subF, 1);                                              % 分 k 类（由objfun决定）
cSearchAgents = floor(SearchAgents / k);    % 并行智能体个数
cFitness = zeros(k, cSearchAgents);               % 保存每类的适应度
cPositions = cell(k, 1);                                     % 保存每类的位置信息
cTracks = cell(k, 1);                                          % 保存每类的航迹信息
cProbPoints = cell(k, 1);                                  % 保存每类的有问题航迹点
cData = cell(k, 1);                                             % 存储每类的检测报告

% 排序
[~, Index] = sort(o_subF, 2, "ascend") ;     % 沿维度2升序排序，返回新矩阵和序号
                                                                       % fitness越小越好
% 聚类
for i = 1:k
    Positions = zeros(cSearchAgents, Dim);
    batchTrack = cell(cSearchAgents, 1);
    batchProbPoints = cell(cSearchAgents, 1);
    batchData = cell(cSearchAgents, 1);
    for j = 1:cSearchAgents
        idx = Index(i, j);
        cFitness(i, j) = o_Fitness(idx); %mat-vector
        Positions(j, :) = WolfPops.Pos(idx, :); %mat-mat
        batchTrack{j} = Tracks{idx}; %cell-cell
        batchProbPoints{j} = o_ProbPoints{idx}; %cell-cell
        batchData{j} = o_Data{idx}; %cell-cell
    end
    cPositions(i) = {Positions}; %cell-mat
    cTracks{i} = batchTrack; %cell-cell
    cProbPoints{i} = batchProbPoints; %cell-cell
    cData{i} = batchData; %cell-cell
end

% 输出
ClassPops.Pos = cPositions;
ClassPops.Tracks = cTracks;
ClassPops.ProbPoints = cProbPoints;
ClassPops.Data = cData;
ClassPops.Fitness = cFitness;
ClassPops.SearchAgents = cSearchAgents;
ClassPops.k = k;
end

repo: Grey-Wolf-Optimizer-for-Path-Planning


text:
function [throttle_factor, counter_out] = manage_clutchgearshift(trigger_gearshift, counter_in, throttle_shiftbreak_s, tSPhysics)
%% Documentation
%
% Author:       Leonhard Hermansdorfer (leo.hermansdorfer@tum.de)
%
% Start Date:   15.09.2021
%
% Description:  this function applies throttle lifting if a gear shift happend.
%               This aims at increasing the vehicle stability during shifting.
%
% Inputs:
%   - trigger_gearshift         gear shift trigger (1 if shift happend, otherwise 0)
%   - counter_in                same value as counter_out from last timestep
%   - throttle_shiftbreak_s     time duration of throttle lift after gear shift
%
% Outputs:
%   - throttle_factor           multiply factor for throttle signal (either 1 or 0)
%   - counter_out               counter which is fed back into the function in the next timestep

    % detect gear shift (direction independent)
    trigger_gearshift = min(1, abs(trigger_gearshift));

    % start counter if a gear shift has happend
    if trigger_gearshift == 1
        counter_in = throttle_shiftbreak_s;
    end

    % count down to zero
    counter_out = counter_in - tSPhysics;

    % calculate throttle factor
    if counter_out > 0
        throttle_factor = 0;
    else
        throttle_factor = 1;
    end

end

author: TUMFTM


text:
function niederreiter2_test03 ( )

%*****************************************************************************80
%
%% NIEDERREITER2_TEST03 tests EXOR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NIEDERREITER2_TEST03\n' );
  fprintf ( 1, '  EXOR returns the exclusive OR of two integers.\n' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '     I     J     EXOR(I,J)\n' );
  fprintf ( 1, '\n' );

  seed = get_seed ( );

  for test = 1 : 10

    [ i, seed ] = i4_uniform ( 0, 100, seed );
    [ j, seed ] = i4_uniform ( 0, 100, seed );
    k = exor ( i, j );

    fprintf ( 1, '%6d %6d %6d\n', i, j, k );

  end

  return
end


text:
%   CLASS Go_Settings
% =========================================================================
%
% DESCRIPTION
%   Class to store all the Processing parameters
%
% EXAMPLE
%   state = Go_Settings();
%
% FOR A LIST OF CONSTANTs and METHODS use doc Go_Settings

%--------------------------------------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:        Andrea Gatti
%  Contributors:      Andrea Gatti, Giulio Taliaferro, ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

classdef Go_Settings < Settings_Interface

    properties (Constant, Access = 'protected')
    end

    % Real constant
    properties(Constant, Access = 'private')
        CUR_INI = 'goGPS_settings.ini';

        LOG_DEFAULT_MODE = 1; % 0 text mode, 1 graphic mode
        LOG_COLOR_MODE = 1;

        GUI_DEFAULT_MODE = 'dark';
        GUI_DEFAULT_EXPORT_MODE = 'light';

        FLAG_EXPORT_TRANSPARENT = true;
    end

    properties (Constant, Access = 'public')
    end

    properties (SetAccess = protected, GetAccess = protected)
        % Location of the current ini file
        log_default_mode = Go_Settings.LOG_DEFAULT_MODE; % 0 text mode, 1 graphic mode
        log_color_mode = Go_Settings.LOG_COLOR_MODE;

        gui_default_mode = Go_Settings.GUI_DEFAULT_MODE;
        gui_default_export_mode = Go_Settings.GUI_DEFAULT_EXPORT_MODE;

        flag_export_transparent = Go_Settings.FLAG_EXPORT_TRANSPARENT;
    end

    properties (SetAccess = public, GetAccess = public)
        cur_ini = Go_Settings.CUR_INI;
        creation_time = GPS_Time(now);
    end

    % =========================================================================
    %%  INIT
    % =========================================================================
    methods(Static,Access = private)
        % Guard the constructor against external invocation.  We only want
        % to allow a single instance of this class.  See description in
        % Singleton superclass.

        function this = Go_Settings(ini_settings_file)
            % Creator
            %
            % SYNTAX
            %   s_obj = Go_Settings(<ini_settings_file>, <prj_home>);

            log = Core.getLogger();
            log.addMarkedMessage('Preparing goGPS settings...');
            log.newLine();

            if nargin == 1 && not(isempty(ini_settings_file))
                this.cur_ini = ini_settings_file;
            end

            if (exist(this.cur_ini, 'file') == 2)
                this.import(this.cur_ini);
            else log.addMarkedMessage('Using default settings');
                log.newLine();
            end
        end
    end

    methods (Static)
        function this = getInstance(ini_settings_file)
            % Concrete implementation.  See Singleton superclass.
            persistent unique_instance_gocfg__
            if isempty(unique_instance_gocfg__)
                if nargin == 1
                    this = Go_Settings(ini_settings_file);
                else
                    this = Go_Settings();
                end
                unique_instance_gocfg__ = this;
            else
                this = unique_instance_gocfg__;
                if nargin == 1
                    this.import(ini_settings_file)
                end
            end
        end
    end

    % =========================================================================
    %%  INTERFACE REQUIREMENTS
    % =========================================================================
    methods (Access = 'public')

        function releoad(this)
            this.import();
        end

        function importIniFile(this, file_name)
            this.import(file_name);
        end

        function import(this, file_name)
            % This function import processing settings from another setting object or ini file
            %
            % SYNTAX
            %   s_obj.import(state)

            if nargin < 2 || isempty(file_name)
                file_name = this.cur_ini;
            end

            ini = Ini_Manager(file_name);

            try this.log_default_mode = ini.getData('LOG', 'default_mode'); catch, keyboard, end
            try this.log_color_mode = ini.getData('LOG', 'color_mode'); catch, keyboard, end

            try this.gui_default_mode = ini.getData('GUI', 'default_mode'); catch, keyboard, end
            try this.gui_default_export_mode = ini.getData('GUI', 'default_export_mode'); catch, keyboard, end

            try this.flag_export_transparent = ini.getData('EXPORT', 'flag_export_transparent'); catch, keyboard, end

            this.check(); % check after import
        end

        function str_cell = export(this, str_cell)
            % Conversion to string ini format of the minimal information needed to reconstruct the this
            %
            % SYNTAX
            %   s_obj.export(str_cell)

            log.addMarkedMessage('Export of goGPS settings - to do...');
        end

        function toString(this)
            ini = Ini_Manager(this.cur_ini);
            ini.showData;
        end
    end

    %%  CHECKING FUNCTIONS - to be kept here to access private parameters
    % =========================================================================
    methods (Access = 'protected')
        function checkLogicalField(this, field_name)
            % Check if a logical field of the object is a valid logical number
            % To make the function works it is needed to have defined the default
            % value of the field as a constant with same name but upper case
            %
            % SYNTAX
            %   this.checkLogicalField(string_field_name);
            this.(field_name) = this.checkLogical(field_name, this.(field_name), this.(upper(field_name)));
        end

        function checkCellStringField(this, field_name, empty_is_valid, check_existence)
            % Check if a string field of the object is a valid string
            % To make the function works it is needed to have defined the default
            % value of the field as a constant with same name but upper case
            %
            % SYNTAX
            %   this.Field(string_field_name, <empty_is_valid == false>, <check_existence == false>);
            switch nargin
                case 2, this.(field_name) = this.checkCellString(field_name, this.(field_name), this.(upper(field_name)));
                case 3, this.(field_name) = this.checkCellString(field_name, this.(field_name), this.(upper(field_name)), empty_is_valid);
                case 4, this.(field_name) = this.checkCellString(field_name, this.(field_name), this.(upper(field_name)), empty_is_valid, check_existence);
                otherwise, error('Settings checkCellStringField called with the wrong number of parameters');
            end
        end

        function is_existing = checkStringField(this, field_name, empty_is_valid, check_existence)
            % Check if a string field of the object is a valid string
            % To make the function works it is needed to have defined the default
            % value of the field as a constant with same name but upper case
            %
            % check_existence 0: do not check
            %                 1: check and do nothing
            %                 2: check and                                try                                                                to correct
            %
            % SYNTAX
            %
            %   this.checkStringField(string_field_name, <empty_is_valid == false>, <check_existence == false>);
            switch nargin
                case 2, [this.(field_name), is_existing] = this.checkString(field_name, this.(field_name), this.(upper(field_name)));
                case 3, [this.(field_name), is_existing] = this.checkString(field_name, this.(field_name), this.(upper(field_name)), empty_is_valid);
                case 4, [this.(field_name), is_existing] = this.checkString(field_name, this.(field_name), this.(upper(field_name)), empty_is_valid, check_existence);
                otherwise, error('Settings checkStringField called with the wrong number of parameters');
            end
        end

        function is_existing = checkPathField(this, field_name, empty_is_valid, check_existence)
            % Check if a string path field of the object is a valid path
            % To make the function works it is needed to have defined the default
            % value of the field as a constant with same name but upper case
            %
            % SYNTAX
            %   this.checkPathField(string_field_name, <empty_is_valid == false>, <check_existence == false>);
            fnp = File_Name_Processor();
            tmp_path = fnp.checkPath(this.(field_name), this.getHomeDir);
            tmp_default_path = fnp.getFullDirPath(fnp.checkPath(this.(upper(field_name)), this.getHomeDir), this.getHomeDir);
            if ~isempty(tmp_path) || empty_is_valid
                this.(field_name) = fnp.getFullDirPath(tmp_path, this.prj_home, [], fnp.getFullDirPath(tmp_default_path));
                switch nargin
                    case 2, [this.(field_name), is_existing] = this.checkString(field_name, tmp_path, tmp_default_path);
                    case 3, [this.(field_name), is_existing] = this.checkString(field_name, tmp_path, tmp_default_path, empty_is_valid);
                    case 4, [this.(field_name), is_existing] = this.checkString(field_name, tmp_path, tmp_default_path, empty_is_valid, check_existence);
                    otherwise, error('Settings checkStringField called with the wrong number of parameters');
                end
                this.(field_name) = fnp.getFullDirPath(this.(field_name), this.getHomeDir);
            else
                this.(field_name) = tmp_default_path;
            end
        end

        function checkNumericField(this, field_name, limits, valid_val)
            % Check if a numeric field of the object is valid
            % To make the function works it is needed to have defined the default
            % value of the field as a constant with same name but upper case
            %
            % SYNTAX
            %   this.checkNumericField(string_field_name, <limits>, <valid_values>);
            switch nargin
                case 2, this.(field_name) = this.checkNumber(field_name, this.(field_name), this.(upper(field_name)));
                case 3, this.(field_name) = this.checkNumber(field_name, this.(field_name), this.(upper(field_name)), limits);
                case 4, this.(field_name) = this.checkNumber(field_name, this.(field_name), this.(upper(field_name)), limits, valid_val);
                otherwise, error('Settings checkNumericField called with the wrong number of parameters');
            end
        end
    end

    % =========================================================================
    %%  TEST PARAMETERS VALIDITY
    % =========================================================================

    methods (Access = 'public')
        function check(this)
            % Check the validity of the fields
            %
            % SYNTAX
            %   this.check();

            this.checkNumericField('log_default_mode', [0 1]);
            this.checkLogicalField('log_color_mode');
            this.checkLogicalField('flag_export_transparent');

            if isdeployed
                this.log_color_mode = false;
            end
        end
    end

    % =========================================================================
    %%  GETTERS
    % =========================================================================
    methods
        function out = getLogMode(this)
            out = this.log_default_mode;
        end

        function out = isLogColorMode(this)
            out = logical(this.log_color_mode);
        end

        function out = getGUIMode(this)
            out = this.gui_default_mode;
        end

        function out = getGUIModeExport(this)
            out = this.gui_default_export_mode;
        end

        function out = isExportTransparent(this)
            out = logical(this.flag_export_transparent);
        end
    end

    %%  SETTERS
    % =========================================================================
    methods
        function out = setLogMode(this, val)
            this.log_default_mode = val;
        end

        function out = setLogColorMode(this, val)
            this.log_color_mode = logical(val);
        end

        function out = setGUIMode(this, val)
            this.gui_default_mode = val;
        end

        function out = setGUIModeExport(this, val)
            this.gui_default_export_mode = val;
        end

        function out = setExportTransparent(this, val)
            this.flag_export_transparent = logical(val);
        end
    end

end

author: goGPS-Project


text:
%function gppca_demo_poster

% Good seeds: 16, 20
seed = 16
randn('state', seed)
rand('state', seed)

inX = 100*(1:200);
inW = 10*rand(2,30);

N = size(inX,2);
M = size(inW,2);

D = 0;

% Trend
D = D + 1;
covfuncX{D} = @gpcov;
logthetaX{D} = log(365*(10+10*rand)); % 10-20 years
initthetaX{D} = logthetaX{D} + 1;%log(365*(10+10*rand));
densityX(D) = 1;

% Periodical (it is difficult to learn from a badly initialized
% periodicity :( )
D = D + 1;
covfuncX{D} = {@gpcovProduct, @gpcov, @gpcovPeriodic};
%logthetaX{D} = log([5*365, 2, 4*365]); % 1 year
logthetaX{D} = log([250*365, 1.6, 3*365]); % decay 100 years, period 3 year
initthetaX{D} = logthetaX{D};%log([300*365, 2, 365*(3.5+1*rand)]);
densityX(D) = 1;

% Short-scale
D = D + 1;
covfuncX{D} = @gpcov;
logthetaX{D} = log(365*1); % 1 year
initthetaX{D} = logthetaX{D} - 0.1;%log(30*(20+20*rand));
densityX(D) = 1;

% Short-scale, almosta noise
D = D + 1;
covfuncX{D} = @gpcovPP;
logthetaX{D} = log(30*4); % 4 months
initthetaX{D} = logthetaX{D} - 0.1;%log(30*(3+6*rand));
densityX(D) = 1;

% Generate latent signals
X = zeros(D,N);
for d=1:D
  X(d,:) = gprnd(inX, logthetaX{d}, covfuncX{d});
end

scales = [2^2 3^2 3^2 1.0^2];
lengthscales = [5 1.5 1.5 0.9];
bias =   [3 0 0 0];
% Working set
% $$$ scales = [2^2 3^2 5^2 1.0^2];
% $$$ lengthscales = [5 2.5 1.5 0.9];
% $$$ bias =   [3 0 0 0];
W = zeros(M,D);
[coord1,coord2] = meshgrid(linspace(0,10,50), linspace(0,10,50));
inWh = [coord1(:)'; coord2(:)'];
W_function = zeros(numel(coord1), D);
for d=1:D
% $$$   covfuncW{d} = {@(logtheta,x1,x2) gpcov(logtheta,x1,x2,@sqdistEuclidean)};
% $$$   logthetaW{d} = [log(2)];
% $$$   initthetaW{d} = [log(6)];
  covfuncW{d} = {@gpcovScale, @(logtheta,x1,x2) gpcov(logtheta,x1,x2,@ ...
                                                    sqdistEuclidean)};
  logthetaW{d} = [log(scales(d)); log(lengthscales(d))];
  initthetaW{d} = logthetaW{d} - 0.1;%[log(5); log(6)];
  densityW(d) = 1;
  W(:,d) = bias(d) + gprnd(inW, logthetaW{d}, covfuncW{d});
% $$$   tmp = bias(d) + gprnd([inW, inWh], logthetaW{d}, covfuncW{d});
% $$$   W(:,d) = tmp(1:M);
% $$$   W_function(:,d) = tmp((M+1):end);
end
for d=1:D
  %tmp = bias(d) + gprnd([inW, inWh], logthetaW{d}, covfuncW{d});
  W_function(:,d) = bias(d) + gppred(inW, W(:,d), 0, inWh, logthetaW{d}, ...
                                     covfuncW{d});
end

% $$$ size_inW = size(inW)
% $$$ size_W = size(W)

% Scale X to unit second moment, so comparison is more valid
scale = sqrt(mean(X.^2,2));
scale = [-1; -1; -1; 1] .* scale; % comparable sign rotation
X = diag(1./scale) * X;
W = W * diag(scale);
W_function = W_function * diag(scale);

%hax = tsplot(inX, X, 'k');
%return

% Data
Y = W*X;
[M,N] = size(Y);

% $$$ variance_noiseless = [(1:M)', std(Y,1,2)]
% $$$ mean_std_noiseless = sqrt(mean(var(Y,1,2)))
% $$$ return

% Noise
Yn = Y + 1*randn(size(Y));

M = size(W,1);
random = randperm(M);
%examples = random([4 9 14 15]);
examples = [19 1 20 5];

% $$$ tsplot(X)
% $$$ tsplot(W')
% $$$ tsplot(Y(examples,:));
% $$$ return

% Missing values
Ynm = Yn;
Itrain = rand(size(Yn))<0.9;
Ynm(Itrain) = nan;
Ytest = Yn;
Ytest(~Itrain) = nan;
Ytest(rand(size(Ytest))<0.87) = nan;

% Create a gap
Ygap = nan*Ynm;
Ynm(:,100:150) = nan;

% $$$ size_Y = size(Ynm)
% $$$ num_of_observations = sum(~isnan(Ynm(:)))
% $$$
% $$$ first_Y = Y(:,1:5)

num_of_observations = sum(sum( ~isnan(Ynm) ))
%return

% $$$ vbpca = false;
% $$$ if vbpca
% $$$   Qvb = pca_full(Ynm,D,'maxiters',50,'rotate2pca',true, 'algorithm','vb');
% $$$
% $$$   % Plot VBPCA latent components
% $$$   varS = zeros(size(Qvb.S));
% $$$   for n=1:N
% $$$     varS(:,n) = diag(Qvb.Sv{n});
% $$$   end
% $$$   hax = tsgpplot(inX', Qvb.S', 2*sqrt(varS)');
% $$$   for i=1:length(hax)
% $$$     set(hax(i), 'xtick', [], 'ytick', []);
% $$$     axes(hax(i));
% $$$     line(100*[100, 100], [-1000 1000], 'Color', [0 0 0])
% $$$     line(100*[150, 150], [-1000 1000], 'Color', [0 0 0])
% $$$     lab = sprintf('x_{%d}(t)', i);
% $$$     ylabel(lab);
% $$$   end
% $$$   xlabel('time, t');
% $$$   set_figure_size(gcf, 12, 7);
% $$$   set_label_fontsize(hax, 6);
% $$$   print('-depsc', '/home/jluttine/thesis_slides/fig_artificial_vbpca_latent');
% $$$
% $$$   % Plot VBPCA predictive distribution
% $$$   Yvb = bsxfun(@plus, Qvb.A*Qvb.S, Qvb.Mu);
% $$$   varYpvb = zeros(size(Yvb));
% $$$   for m=1:M
% $$$     CovA = Qvb.Av{m};
% $$$     for n=1:N
% $$$       CovS = Qvb.Sv{n};
% $$$       varYvb(m,n) = Qvb.A(m,:)*CovS*Qvb.A(m,:)' + Qvb.S(:,n)'*CovA*Qvb.S(:,n) ...
% $$$           + traceprod(CovS,CovA,true);
% $$$     end
% $$$   end
% $$$   varYvb = bsxfun(@plus, varYvb, Qvb.Muv + Qvb.V);
% $$$   hax = tsgpplot(inX', Yvb(examples,:)', 2*sqrt(varYvb(examples,:))');
% $$$   addtsplot(inX, Yn(examples,:), 'r-');
% $$$   addtsplot(inX, Yvb(examples,:), 'k-'); % draw vbpca mean again..
% $$$   addtsplot(inX, Ynm(examples,:), '+', 'MarkerSize', 7, 'Color', [0 0 1]);
% $$$   for i=1:length(hax)
% $$$     set(hax(i), 'xtick', [], 'ytick', []);
% $$$     axes(hax(i));
% $$$     line(100*[100, 100], [-1000 1000], 'Color', [0 0 0])
% $$$     line(100*[150, 150], [-1000 1000], 'Color', [0 0 0])
% $$$     lab = sprintf('y_{%d}(t)', examples(i));
% $$$     ylabel(lab);
% $$$   end
% $$$   %xlabel('time, t');
% $$$   set_figure_size(gcf, 12, 10.5);
% $$$   set_label_fontsize(hax, 6);
% $$$   print('-depsc', ['/home/jluttine/papers/2009NIPS/poster/' ...
% $$$                     'fig_artificial_vbpca_predictive']);
% $$$ end

%return

Qgp = vbgppcamv_full(Ynm, D,...
                     inW, inX,...
                     covfuncW, initthetaW,...
                     covfuncX,initthetaX, ...
                     'maxiter',50, ...
                     'pseudodensityx', densityX, ...
                     'pseudodensityw', densityW, ...
                     'loglikelihood', true, ...
                     'updatehyper', [5 10 20 50 80 100], ...
                     'updatepseudox', false, ...
                     'updatepseudow', false, ...
                     'maxsearchx', 3, ...
                     'maxsearchw', 3, ...
                     'checkgradw', false, ...
                     'checkgradx', false);

Ygp = Qgp.W * Qgp.X;
if false
  disp('Using inaccurate prediction variance.');
  varYgp = Qgp.W.^2*Qgp.varX + Qgp.varW*Qgp.X.^2 + Qgp.varW*Qgp.varX + ...
           1/Qgp.tau;
else
  disp('Using accurate prediction variance.')
  varYgp = zeros(size(Ygp));
  for m=1:M
    CovWm = Qgp.CovW(Qgp.indsW(m,:),Qgp.indsW(m,:));
    for n=1:N
      CovXn = Qgp.CovX(Qgp.indsX(:,n),Qgp.indsX(:,n));
      varYgp(m,n) = Qgp.W(m,:)*CovXn*Qgp.W(m,:)' + Qgp.X(:,n)'*CovWm*Qgp.X(:,n) ...
          + traceprod(CovXn,CovWm,true);
    end
  end
  varYgp = varYgp + 1/Qgp.tau;
end

% Plot latent signals
hax = tsgpplot(inX', Qgp.X', 2*sqrt(Qgp.varX)');
for i=1:length(hax)
  set(hax(i), 'xtick', [], 'ytick', []);
  axes(hax(i));
  line(100*[100, 100], [-1000 1000], 'Color', [0 0 0])
  line(100*[150, 150], [-1000 1000], 'Color', [0 0 0])
  lab = sprintf('x_{%d}(t)', i);
%  ylabel(lab);
end
%xlabel('time, t');
set_subplot_positions(hax, 4, 1, [0.01 0.01 0.01 0.01], [0.02 0.02]);
set_figure_size(gcf, 7, 6);
set_label_fontsize(hax, 6);
set(gcf, 'Color','none');
set(hax, 'Color', 'none');
export_fig('/home/jluttine/thesis/slides/novac2010_artificial_latent', ...
           '-eps');
% $$$ print('-depsc', '/home/jluttine/papers/2009NIPS/poster/fig_artificial_latent');

% $$$ % Plot true latent signals
% $$$ hax = tsplot(inX, X, 'k');
% $$$ hold on
% $$$ yl = max(-min(X,[],2), max(X,[],2))
% $$$ for i=1:length(hax)
% $$$   set(hax(i), 'YLim', [-yl(i) yl(i)]);
% $$$   set(hax(i), 'xtick', [], 'ytick', []);
% $$$   axes(hax(i));
% $$$   line(100*[100, 100], [-1000 1000], 'Color', [0 0 0])
% $$$   line(100*[150, 150], [-1000 1000], 'Color', [0 0 0])
% $$$   lab = sprintf('x_{%d}(t)', i);
% $$$   ylabel(lab);
% $$$ end
% $$$ xlabel('time, t');
% $$$ set_figure_size(gcf, 7, 6);
% $$$ set_label_fontsize(hax, 6);
% $$$ print('-depsc', '/home/jluttine/papers/2009NIPS/poster/fig_artificial_true_latent');

% Compare latent signals
% $$$ Qgp.X = bsxfun(@minus, Qgp.X, mean(Qgp.X,2));
% $$$ Qgp.varX = bsxfun(@rdivide, Qgp.varX, std(Qgp.X,1,2).^2);
% $$$ Qgp.X = bsxfun(@rdivide, Qgp.X, std(Qgp.X,1,2));
% $$$ X = bsxfun(@minus, X, mean(X,2));
% $$$ X = bsxfun(@rdivide, X, std(X,1,2));
% $$$ tsgpplot(inX', Qgp.X', 2*sqrt(Qgp.varX)');
% $$$ addtsplot(inX, X, 'r')

% Show predictive distribution
hax = tsgpplot(inX', Ygp(examples,:)', 2*sqrt(varYgp(examples,:))');
addtsplot(inX, Yn(examples,:), 'r-');%, 'MarkerSize', 1)
addtsplot(inX, Ygp(examples,:), 'k-') % draw gp mean again..
addtsplot(inX, Ynm(examples,:), '+', 'MarkerSize', 4, 'Color', [0 0 1]);
% $$$ addtsplot(inX, Ytest(examples,:), 'o', 'MarkerSize', 4, 'Color', [0 ...
% $$$                     0 0.8]);
for i=1:length(hax)
  set(hax(i), 'xtick', [], 'ytick', []);
  axes(hax(i));
  line(100*[100, 100], [-1000 1000], 'Color', [0 0 0])
  line(100*[150, 150], [-1000 1000], 'Color', [0 0 0])
% $$$   lab = sprintf('y_{%d}(t)', examples(i));
% $$$   ylabel(lab);
end
%xlabel('time, t');
set_subplot_positions(hax, 4, 1, [0.01 0.01 0.01 0.01], [0.02 0.02]);
set_figure_size(gcf, 7, 6);
set_label_fontsize(hax, 6);
set(gcf, 'Color','none');
set(hax, 'Color', 'none');
export_fig('/home/jluttine/thesis/slides/novac2010_artificial_predictive', ...
           '-eps');
% $$$ print('-depsc', '/home/jluttine/papers/2009NIPS/poster/fig_artificial_predictive');

exp(Qgp.logthetaX{1})
exp(Qgp.logthetaX{2})
exp(Qgp.logthetaX{3})

%return

% Plot spatial components
Mh = size(inWh,2);
Wh = zeros(Mh,D);
varWh = zeros(Mh,D);
fig1 = figure;
% $$$ fig2 = figure;
% $$$ fig3 = figure;
hax = [];
for d=1:D
  %first = (d-1)*M + 1;
  %last = first + M - 1;
  %inds = first:last;
  inds = Qgp.indsW(:,d);
  [Wh(:,d), varWh(:,d)] = gppred(Qgp.inW, Qgp.W(:,d), Qgp.CovW(inds,inds), ...
                                 inWh, Qgp.logthetaW{d}, Qgp.covfuncW{d});

  % Plot mean map
  figure(fig1)
  subplot(2,ceil(D/2),d)
  contourf(coord1,coord2,reshape(Wh(:,d), size(coord1)),20);
  hax(d,1) = gca;
  hold on
  plot(inW(1,:),inW(2,:), 'kx', 'MarkerSize', 8);
  plot(inW(1,examples),inW(2,examples), 'ko', 'MarkerSize', 8);
  %  plot(inW(1,examples),inW(2,examples), 'kx', 'MarkerSize', 10);
  shading('flat');
  cl = get(gca, 'clim');
  cl = max(abs(cl));
  set(gca, 'clim', [-cl cl]);
  mapcolormap;
  set(gca, 'xtick', [], 'ytick', []);
  %pbaspect([1 1 1]);
%  h_cb(d,1) = colorbar('SouthOutside');

% $$$   % Plot uncertainty map
% $$$   figure(fig2)
% $$$ %  subplot(1,D,d);
% $$$   subplot(2,D,D+d);
% $$$   contourf(coord1,coord2,reshape(sqrt(varWh(:,d)), size(coord1)),20);
% $$$   hax(d,2) = gca;
% $$$   hold on
% $$$   plot(inW(1,:),inW(2,:), 'kx', 'MarkerSize', 8);
% $$$   plot(inW(1,examples),inW(2,examples), 'ko', 'MarkerSize', 8);
% $$$   %  plot(inW(1,examples),inW(2,examples), 'kx', 'MarkerSize', 10);
% $$$   shading('flat');
% $$$   cl = get(gca, 'clim');
% $$$   cl = max(abs(cl));
% $$$   set(gca, 'clim', [0 cl]);
% $$$   %mapcolormap;
% $$$   m = colormap('hot');
% $$$   colormap(m(end:-1:1,:));
% $$$   set(gca, 'xtick', [], 'ytick', []);
% $$$   %pbaspect([1 1 1]);
% $$$   h_cb(d,2) = colorbar('SouthOutside');

% $$$   % Plot real map
% $$$   figure(fig3)
% $$$   subplot(1,D,d)
% $$$   contourf(coord1,coord2,reshape(W_function(:,d), size(coord1)),20);
% $$$   hax(d,3) = gca;
% $$$   hold on
% $$$   plot(inW(1,:),inW(2,:), 'kx', 'MarkerSize', 8);
% $$$   plot(inW(1,examples),inW(2,examples), 'ko', 'MarkerSize', 8);
% $$$   %  plot(inW(1,examples),inW(2,examples), 'kx', 'MarkerSize', 10);
% $$$   shading('flat');
% $$$   cl = get(gca, 'clim');
% $$$   cl = max(abs(cl));
% $$$   set(gca, 'clim', [-cl cl]);
% $$$   mapcolormap;
% $$$   set(gca, 'xtick', [], 'ytick', []);
% $$$   %pbaspect([1 1 1]);
% $$$   h_cb(d,3) = colorbar('SouthOutside');

end
cb_height = 0.0;
cb_textheight = 0.0;
cb_sep = 0.0;

for i=1:1
  set_subplot_positions(hax(:,i), 2, 2, [0.01 0.01 0.01 0.01], [0.02 0.02]);
end
for d=1:numel(hax)
  pos_ax = get( hax(d), 'Position');
  cb_pos = pos_ax;
  cb_pos(2) = pos_ax(2) - cb_height - cb_sep;
  cb_pos(4) = cb_height;
  axis(hax(d), 'square')
%  set(h_cb(d), 'Position', cb_pos, 'FontSize', 7);
end

set_figure_size(fig1, 6, 6);
% $$$ set_figure_size(fig2, 12, 4);
% $$$ set_figure_size(fig3, 12, 4);
print(fig1, '-depsc', ['/home/jluttine/thesis/slides/' ...
                    'novac2010_artificial_loadings']);
% $$$ print(fig2, '-depsc', ['/home/jluttine/papers/2009NIPS/poster/' ...
% $$$                     'novac2010_artificial_loadings_uncertainty']);
% $$$ print(fig3, '-depsc', ['/home/jluttine/papers/2009NIPS/poster/' ...
% $$$                     'novac2010_artificial_true_loadings']);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%
% 	System response

% Conversion between continuous time and discrete time systems

%continious to discrete
n=[1 10];
d=[1 0];
Hs=tf(n,d)

Ts=0.2;
Hz=c2d(Hs,Ts,'zoh')

%discrete to continious
num= [1 1];
den=[1 -1];
Ts=0.2;
Hz=tf(num,den,Ts)

Hs=d2c(Hz,'zoh')


text:
function selected_indices=cosmo_anova_feature_selector(dataset, how_many)
% find the features that show the most variance between classes
%
% selected_indices=cosmo_anova_feature_selector(dataset, how_many)
%
% Inputs:
%  dataset          struct with .samples and .sa.targets
%  how_many         value between 0 and 1 keeps how_many*100% features;
%                   values >=1 keeps how_many features
%
% Output:
%  selected_indices   feature ids in dataset with most variance between
%                     classes.
%
% Example:
%     ds=cosmo_synthetic_dataset();
%     disp(size(ds.samples))
%     %|| [ 6 6 ]
%     cosmo_anova_feature_selector(ds,.45) % find best ~45% of features
%     %|| [ 2 4 5 ]
%     cosmo_anova_feature_selector(ds,4) % find best 4 features
%     %|| [ 2 4 5 3 ]
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    fstat=cosmo_stat(dataset,'F');
    fvalues=fstat.samples;

    % ensure that nan values are not selected by setting them to
    % an impossible low F value
    fvalues(isnan(fvalues))=-1;

    % sort by F values, largest first
    [unused, idxs]=sort(fvalues,'descend');

    % determine features to select
    nfeatures=size(dataset.samples,2);

    if how_many>=1
        if round(how_many)~=how_many
            error('how_many>=1 is not an integer');
        elseif how_many>nfeatures
            error('dataset has %d features, cannot return %d',...
                    nfeatures,how_many);
        end
        nkeep=how_many;
    else
        nkeep=round(how_many*nfeatures);
    end

    selected_indices=idxs(1:nkeep);

    % throw an error if any indices with NaN F values
    if any(fvalues(selected_indices)<0)
        idx=find(fvalues(selected_indices)<0,1);
        error('Feature %d has NaN Fscore', selected_indices(idx));
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni

%problem 2 - Fourier Transform of sinc(t)

syms t w
x=sin(pi*t)/(pi*t);
X=fourier(x,w);
ezplot(X, [-10 10])
legend('X(\Omega)');
xlabel('\Omega');


text:
close all
clear all
image = im2double(imread('..\2.tif'));

%% proof that txt files and stage 4 size is not equal.
prophoto2 = im2double(imread('..\3.tif'));
%center aligned
eucliudian_error1 = sqrt((image(:,:,1)-...
    prophoto2(:,:,1)).^2 + ...
    (image(:,:,2)-...
    prophoto2(:,:,2)).^2 + ...
    (image(:,:,3)-...
    prophoto2(:,:,3)).^2);
sum(eucliudian_error1(:))
figure,
imagesc(eucliudian_error1, [0 0.01]);

% abs()

author: karaimer

repo: camera-pipeline-UI


text:
clear all; close all; clc;

close all
figure(1)
load catData_w.mat
load dogData_w.mat
CD=[dog_wave cat_wave];
[u,s,v]=svd(CD-mean(CD(:)));

subplot(2,2,1)
plot(v(1:80,2),v(1:80,4),'ro','Linewidth',[1],'MarkerEdgeColor','k',...
                'MarkerFaceColor',[0 1 0.2],...
                'MarkerSize',8), hold on   % [0.49 1 .63]
plot(v(81:end,2),v(81:end,4),'bo','Linewidth',[1],'MarkerEdgeColor','k',...
                'MarkerFaceColor',[0.9 0 1],...
                'MarkerSize',8)
set(gca,'Fontsize',[15]), hold on

dogcat=v(:,2:2:4);
GMModel=fitgmdist(dogcat,2)
AIC= GMModel.AIC

subplot(2,2,1)
h=ezcontour(@(x1,x2)pdf(GMModel,[x1 x2]));
xlabel(''),ylabel(''),title('')
set(h,'Linecolor','k','Linewidth',[2])
axis([-.25 .25 -.25 .25])

subplot(2,2,2)
h=ezmesh(@(x1,x2)pdf(GMModel,[x1 x2]));
xlabel(''),ylabel(''),title(''), view(-16,40)
set(gca,'Fontsize',[15]), colormap(gray)
axis([-.25 .25 -.25 .25 0 30])

%% AIC scores

figure(2)
AIC = zeros(1,4);
GMModels = cell(1,4);
options = statset('MaxIter',500);
for k = 1:4
    GMModels = fitgmdist(dogcat,k,'Options',options,'CovarianceType','diagonal');
    AIC(k)= GMModels.AIC
    subplot(2,2,k)
    h=ezmeshc(@(x1,x2)pdf(GMModels,[x1 x2]));
    xlabel(''),ylabel(''),title('')
%    set(h,'Linecolor','k','Linewidth',[2])
%    axis([-.25 .25 -.25 .25 0 30])
end

author: dynamicslab


text:
function r8sto_print_some ( n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% R8STO_PRINT_SOME prints some of a R8STO matrix.
%
%  Discussion:
%
%    The R8STO storage format is used for a symmetric Toeplitz matrix.
%    It stores the N elements of the first row.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, real A(N), the R8STO matrix.
%
%    Input, integer ILO, JLO, IHI, JHI, the first row and
%    column, and the last row and column to be printed.
%
%    Input, string TITLE, a title.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );

  incx = 5;
%
%  Print the columns of the matrix, in strips of 5.
%
  for j2lo = jlo : incx : jhi

    j2hi = j2lo + incx - 1;
    j2hi = min ( j2hi, n );
    j2hi = min ( j2hi, jhi );

    inc = j2hi + 1 - j2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col: ' );

    for j = j2lo : j2hi
      fprintf ( 1, '%7d       ', j );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row\n' );
    fprintf ( 1, '  ---\n' );
%
%  Determine the range of the rows in this strip.
%
    i2lo = max ( ilo, 1 );
    i2hi = min ( ihi, n );

    for i = i2lo : i2hi

      fprintf ( 1, '%4d  ', i );
%
%  Print out (up to) 5 entries in row I, that lie in the current strip.
%
      for j2 = 1 : inc

        j = j2lo - 1 + j2;

        if ( i <= j )
          aij = a(1+j-i);
        else
          aij = a(1+i-j);
        end

        fprintf ( 1, '%12g  ', aij );

      end

      fprintf ( 1, '\n' );

    end

  end

  return
end


text:
function results = spm_tests(varargin)
% Unit Testing Framework
% FORMAT results = spm_tests(name,value,...)
% name,value  - pairs of optional parameter names and values:
%     verbose:   verbosity level of test run progress report [default: 2]
%     display:   display test results [default: false]
%     coverage:  display code coverage [default: false]
%     cobertura: save code coverage results in the Cobertura XML format [default: false]
%     tag:       test tag selector [default: '', ie all tests]
%     tap:       save a Test Anything Protocol (TAP) file [default: false]
%     test:      name of function to test [default: '', ie all tests]
%
% results     - TestResult array containing information describing the
%               result of running the test suite.
%__________________________________________________________________________
% Copyright (C) 2015-2019 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_tests.m 7548 2019-03-20 10:08:49Z guillaume $

if spm_check_version('matlab','8.3') < 0
    error('Unit Tests require MATLAB R2014a or above.');
end

SVNid = '$Rev: 7548 $';
SPMid = spm('FnBanner',mfilename,SVNid);

%-Input parameters
%--------------------------------------------------------------------------
options = struct('verbose',2, 'display',false, 'coverage',false, ...
                 'cobertura',false, 'tag', '', 'tap',false, 'test','');
if nargin
    if isstruct(varargin{1})
        fn = fieldnames(varargin{1});
        for i=1:numel(fn)
            options.(fn{i}) = varargin{1}.(fn{i});
        end
    else
        for i=1:2:numel(varargin)
            options.(varargin{i}) = varargin{i+1};
        end
    end
end

%-Unit tests
%==========================================================================

%-Create a TestSuite
%--------------------------------------------------------------------------
import matlab.unittest.TestSuite;
import matlab.unittest.selectors.*;
tests = fullfile(spm('Dir'),'tests');
if isempty(options.test)
    suite = TestSuite.fromFolder(tests, 'IncludingSubfolders', true);
else
    suite = [];
    options.test = cellstr(options.test);
    for i=1:numel(options.test)
        prefix = 'test_';
        if strncmp(options.test{i},prefix,length(prefix))
            prefix = '';
        end
        mtest = fullfile(tests,[prefix spm_file(options.test{i},'ext','.m')]);
        if ~spm_existfile(mtest)
            warning('SPM:tests:fileNotFound',...
                'No tests found for %s',options.test{i});
            continue
        end
        suite = [suite, TestSuite.fromFile(mtest)];
        % if i==1
        %     suite = TestSuite.fromFile(mtest);
        % else
        %     suite(i) = TestSuite.fromFile(mtest);
        % end
    end
end
if ~isempty(options.tag)
    suite = suite.selectIf(~HasTag | HasTag(options.tag));
end

%-Create a TestRunner
%--------------------------------------------------------------------------
import matlab.unittest.TestRunner;
import matlab.unittest.plugins.*
if ~options.verbose
    runner = TestRunner.withNoPlugins;
else
    runner = TestRunner.withTextOutput('Verbosity',options.verbose);
end

if options.coverage
    plugin = CodeCoveragePlugin.forFolder(spm('Dir'));
    runner.addPlugin(plugin);
end

if options.cobertura
    d = getenv('WORKSPACE');
    if isempty(d), d = spm('Dir'); end
    coberturaFile = fullfile(d,'spm_CoverageResults.xml');
    plugin = CodeCoveragePlugin.forFolder(spm('Dir'),...
        'Producing',codecoverage.CoberturaFormat(coberturaFile));
    runner.addPlugin(plugin);
end

if options.tap
    d = getenv('WORKSPACE');
    if isempty(d), d = spm('Dir'); end
    tapFile = fullfile(d,'spm_tests.tap');
    plugin = TAPPlugin.producingVersion13(ToFile(tapFile));
    runner.addPlugin(plugin);
end

%-Run tests
%--------------------------------------------------------------------------
if ~isempty(suite)
    results = runner.run(suite);
else
    results = struct('Passed',{},'Failed',{},'Incomplete',{},'Duration',{});
end

%-Display test results
%--------------------------------------------------------------------------
if options.display
    fprintf(['Totals (%d tests):\n\t%d Passed, %d Failed, %d Incomplete.\n' ...
        '\t%f seconds testing time.\n\n'],numel(results),nnz([results.Passed]),...
        nnz([results.Failed]),nnz([results.Incomplete]),sum([results.Duration]));

    disp(table(results));
end

author: spm


text:
function [E, Emin, Emax, Es] = grbm_energy (x, W, vbias, hbias, sigmas)

n_samples = size(x,1);

Wh = ones(n_samples,1)*hbias' + bsxfun(@rdivide, x, sigmas.^2) * W;

Es = bsxfun(@minus,x,vbias');
Es = sum(bsxfun(@rdivide, Es.^2, sigmas.^2),2)/2;

Wh_plus = max(Wh, 0);
Es = Es - sum(log(exp(-Wh_plus)+exp(Wh-Wh_plus))+Wh_plus,2);
Es = Es';

Emin = min(Es);
Emax = max(Es);
E = mean(Es);

end

author: kyunghyuncho

repo: deepmat


text:
% TD | CP-ALS | PARAFAC/CP decomposition solved by Alternating Least Squares
% process_video('TD', 'CP-ALS', 'dataset/demo.avi', 'output/demo_CP-ALS.avi');

r = 10;
A = double(T);
L = double(cp_als(T,r,'dimorder',[3 2 1]));
S = A - L;


text:
% ty_axes stores all information necessary to generate a nice looking
% ty-Plot for possible multiple timeseries datasets. It provides methods
% to generate this plot if a target axes is given.
%
% It expects the following interface:
% ty_axes(ts_list, yaxis_string)
% with ts_list being a cell array of multiple time series to display and
% yaxis_string being the label placed at the y axis.
classdef ty_axes
   properties
    % cell array which holds multiple timeseries, arrangement does not matter
    ts_list
    % stores the y axis label shared by all timeseries
    yaxis_string
   end
   methods
    % constructor
    function obj = ty_axes(ts_list, yaxis_string)
      obj.ts_list = ts_list;
      obj.yaxis_string = yaxis_string;
    end
    % plots all timeseries objects and setups the axes properties right
    function plot_all(obj, ax, tStart, tEnd)
      % setup general properties
      grid(ax, 'on');
      hold(ax, 'on');
      xlabel(ax, 'Time in s');
      ylabel(ax, obj.yaxis_string{1});
      % initialize plot handles which store the resulting plots to match
      % with the legends later
      ph = [];
      legend_bool = [];
      ph_idx = 1;
      % iterate via all timeseries
      for i = 1:1:numel(obj.ts_list)
        if(nargin==4)
          % if time points are given, specify them
          ph = [ph, obj.plot_ts(obj.ts_list{i}, ax, tStart, tEnd)];
          legend_bool(ph_idx) = 1;
        else
          % otherwise plot complete timeseries by suppression of time points
          ph = [ph, obj.plot_ts(obj.ts_list{i}, ax)];
          legend_bool(ph_idx) = 1;
        end
        ph_idx = ph_idx + 1;
      end
      % tidy up
      if(numel(obj.ts_list) > 1)
        % only display legend if more than one signal is present
        legend(ax, ph(logical(legend_bool)), 'Location', 'best');
      end
      box(ax, 'on');
      hold(ax, 'off');
    end
   % helper function which takes care of time management
   function ph = plot_ts(obj, ts, ax, tStart, tEnd)
    if(nargin~=5)
      % if not all aguments are given, assume that the timeseries is plotted completly
      tStart = ts.Time(1);
      tEnd = ts.Time(end);
    end
    try
      [idxStart, idxEnd] = find_ts_idx(ts, tStart, tEnd);
      % plot data and configure axes
      ph = stairs(ax, ts.Time(idxStart:idxEnd), ts.Data(idxStart:idxEnd), 'DisplayName', strrep(ts.Name,'_',' '));
    catch
      warning(['Error during plotting of ' ts.Name]);
    end
   end
 end
end


text:
function ng = ellipsoid_grid_count ( n, r, c )

%*****************************************************************************80
%
%% ELLIPSOID_GRID_COUNT counts the grid points inside an ellipsoid.
%
%  Discussion:
%
%    The ellipsoid is specified as
%
%      ( ( X - C1 ) / R1 )^2
%    + ( ( Y - C2 ) / R2 )^2
%    + ( ( Z - C3 ) / R3 )^2 = 1
%
%    The user supplies a number N.  There will be N+1 grid points along
%    the shortest axis.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of subintervals.
%
%    Input, real R(3), the half axis lengths.
%
%    Input, real C(3), the center of the ellipsoid.
%
%    Output, integer NG, the number of grid points inside the ellipsoid.
%
  if ( r(1) == min ( r ) )
    h = 2 * r(1) / ( 2 * n + 1 );
    ni = n;
    nj = ceil ( r(2) / r(1) ) * n;
    nk = ceil ( r(3) / r(1) ) * n;
  elseif ( r(2) == min ( r ) )
    h = 2 * r(2) / ( 2 * n + 1 );
    nj = n;
    ni = ceil ( r(1) / r(2) ) * n;
    nk = ceil ( r(3) / r(2) ) * n;
  else
    h = 2 * r(3) / ( 2 * n + 1 );
    nk = n;
    ni = ceil ( r(1) / r(3) ) * n;
    nj = ceil ( r(2) / r(3) ) * n;
  end

  ng = 0;

  for k = 0 : nk
    z = c(3) + k * h;
    for j = 0 : nj
      y = c(2) + j * h;
      for i = 0 : ni
        x = c(1) + i * h;
        if ( 1 < ( ( x - c(1) ) / r(1) ).^2 ...
               + ( ( y - c(2) ) / r(2) ).^2 ...
               + ( ( z - c(3) ) / r(3) ).^2 )
          break
        end
%
%  At least one point is generated, but more possible by symmetry.
%
        np = 1;
        if ( 0 < i )
          np = 2 * np;
        end
        if ( 0 < j )
          np = 2 * np;
        end
        if ( 0 < k )
          np = 2 * np;
        end
        ng = ng + np;

      end
    end
  end

  return
end


text:
% Defines a note
%
% Author: Jonathan Karr
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 5/7/2009
classdef Note < edu.stanford.covert.cell.kb.KnowledgeBaseObject
    properties
    end

    methods
        function this = Note(knowledgeBase, wid, wholeCellModelID, name, ...
                comments, crossReferences)

            if nargin == 0; return; end;

            this = edu.stanford.covert.cell.kb.Note.empty(size(wid,1),0);
            this(size(wid,1),1) = edu.stanford.covert.cell.kb.Note;
            for i = 1:size(wid,1)
                this(i,1).idx = i;
                this(i,1).knowledgeBase = knowledgeBase;
                this(i,1).wid = wid(i);
                this(i,1).wholeCellModelID = wholeCellModelID{i};
                if ~isempty(name); this(i,1).name = name{i}; end;
                if exist('comments','var') && ~isempty(comments); this(i,1).comments = comments{i}; end;
                if exist('crossReferences','var')
                    if size(crossReferences,1)>1
                        this(i,1).crossReferences = crossReferences(i);
                    else
                        this(i,1).crossReferences = struct;
                        fields = fieldnames(crossReferences);
                        for j = 1:size(fields,1)
                            values = crossReferences.(fields{j});
                            this(i,1).crossReferences.(fields{j}) = values(i);
                        end
                    end
                end
            end
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [sol, info] = gsp_proj_filterbank(x, ~ , G, W, y, param)
%GSP_PROJ_FILTERBANK Projection onto the synthesis coefficients
%   Usage:  sol = gsp_proj_filterbank(x, 0, G, W, y, param);
%           sol = gsp_proj_filterbank(x, 0, G, W, y);
%           [sol, info] = gsp_proj_filterbank(...)
%
%   Input parameters:
%         x     : Input signal
%         G     : Graph structure
%         W     : Filterbank (cell array of functions)
%         y     : Measurements
%         param : Structure of optional parameters
%   Output parameters
%         sol   : Solution.
%         info  : Structure summarizing informations at convergence
%
%   `gsp_proj_filterbank(x, gamma, G, W, param)` can solves:
%
%   .. sol = argmin_{z} 0.5*||x - z||_2^2  such that W^* x =y
%
%   .. math::  sol = \min_{z} \frac{1}{2} \|x - z\|_2^2 \text{ s. t. }  W^* x = y
%
%   Where $W$ is the linear analysis operator associated with the
%   filterbank.
%
%   The function can use different techniques
%
%   * 'exact' : if the Fourier basis is computed, go for this one
%   * 'cheby' : use the pseudo-inverse filters of the filterbank with
%     chebyshev approximation. It works well for well-conditionned
%     filterbanks.
%   * 'lanczos' : use the pseudo-inverse filters of the filterbank with
%     lanczos approximation. It works well for well-conditionned
%     filterbanks.
%   * 'proj_b2': scallable and robust way to do it. However, the
%     convergence maybe slow and might require a lot of filtering
%     operations.
%
%   param is a Matlab structure containing the following fields:
%
%   * *param.verbose* : 0 no log, 1 a summary at convergence, 2 print main
%     steps (default: 1)
%   * *param.eps* : tolerance for the pseudo inverse method
%   * *param.proj_method*: selected method
%
%   info is a Matlab structure containing the following fields:
%
%   * *info.algo* : Algorithm used
%   * *info.iter* : Number of iteration
%   * *info.time* : Time of exectution of the function in sec.
%   * *info.final_eval* : Final evaluation of the function
%   * *info.crit* : Stopping critterion used
%
%
%   See also:  gsp_solve_l1 gsp_proj_b2_filterbank
%

% Author: Nathanael Perraudin
% Date: 25 March 2014
% Testing: test_gsp_proj_filerbank

if nargin < 5
    error('GSP_PROJ_FILTERBANK: Not enought input arguments');
end

if nargin < 6, param=struct; end

if ~isfield(param, 'eps'), param.eps = 1e-8; end
if ~isfield(param, 'proj_method'),
    if gsp_check_fourier(G)
        param.proj_method = 'exact';
    else
        [A,B] = gsp_filterbank_bounds(G,W);
        if B/A > 20
            warning('Filterbank ill-conditioned, going for a slow method');
            param.proj_method = 'primal_dual';
        else
            param.proj_method = 'cheby';
        end
    end
end

t = tic;

if size(y,2)==1 && size(x,2)>1
    y = repmat(y,1,size(x,2));
end
if isfield(param,'method')
    param2 = rmfield(param,'method');
else
    param2 = param;
end
if strcmp(param.proj_method, 'proj_b2') || strcmp(param.proj_method, 'primal_dual')
    [~,B] = gsp_filterbank_bounds(G,W);
    if ~isfield(param,'paramproj'), param.paramproj = struct; end
    paramproj = param.paramproj;
    paramproj.A = @(x) gsp_filter_synthesis(G,W,x,param2);
    paramproj.At = @(x) gsp_filter_analysis(G,W,x,param2);
    paramproj.nu = B^2;
    paramproj.epsilon = sqrt(G.N)*param.eps;
    paramproj.y = y;
    paramproj.method  = param.proj_method;

    [sol,info] = proj_linear_eq(x,0,paramproj);

else
    if strcmp(param.proj_method, 'cheby')  || strcmp(param.proj_method, 'lanczos')
        [A,B] = gsp_filterbank_bounds(G,W);
        if B/A > 20
            warning('Filterbank ill-conditioned, check your solution!');
        end
    end
    Wd = gsp_design_can_dual(W,param.eps);
    sol =  x - gsp_filter_analysis(G,Wd,(gsp_filter_synthesis(G,W,x,param2)-y),param2);

    info.iter = 1;
    info.final_eval = 0;
    info.crit = 'Direct computation';
end

    info.time = toc(t);
    info.algo = param.proj_method;

end

repo: gspbox


text:
function predicted=cosmo_classify_nn(samples_train, targets_train, samples_test, unused)
% nearest neighbor classifier
%
% predicted=cosmo_classify_nn(samples_train, targets_train, samples_test[, opt])
%
% Inputs:
%   samples_train      PxR training data for P samples and R features
%   targets_train      Px1 training data classes
%   samples_test       QxR test data
%   opt                (currently ignored)
%
% Output:
%   predicted          Qx1 predicted data classes for samples_test
%
% Example:
%     ds=cosmo_synthetic_dataset('ntargets',5,'nchunks',15);
%     test_chunk=1;
%     te=cosmo_slice(ds,ds.sa.chunks==test_chunk);
%     tr=cosmo_slice(ds,ds.sa.chunks~=test_chunk);
%     pred=cosmo_classify_nn(tr.samples,tr.sa.targets,te.samples,struct);
%     % show targets and predicted labels (100% accuracy)
%     disp([te.sa.targets pred])
%     %||      1     1
%     %||      2     2
%     %||      3     3
%     %||      4     4
%     %||      5     5
%
% See also: cosmo_crossvalidate, cosmo_crossvalidation_measure
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    [ntrain, nfeatures]=size(samples_train);
    [ntest, nfeatures_]=size(samples_test);
    ntrain_=numel(targets_train);

    if nfeatures~=nfeatures_ || ntrain_~=ntrain, error('illegal input size'); end

    % allocate space for output
    predicted=zeros(ntest,1);

    for k=1:ntest
        % for each sample in the test set:
        %
        % - compute its squared euclidian distance to each sample in
        %   the train set, and store this in a vector
        %   squared_distances (which must have size ntrain x 1).
        %   For two vectors a=[a_1, a_2, ..., a_N] and b=[b_1, b_2, ..., b_N],
        %   the squared euclidean distance between a and b is:
        %       (a_1 - b_1)^2 + (a_2 - b_2)^2 + ... + (a_N - b_N)^2
        %
        % - assign the class label of the sample in the training set that has
        %   the smallest squared distance.
        %
        % >@@>
        % compute difference to each sample in the training set
        delta=bsxfun(@minus, samples_train, samples_test(k,:));

        % compute distance (sqrt is unnecessary because monotonic)
        squared_distances=sum(delta.^2,2);

        % the following code is equivalent to (but slower than) the code above:
        %   squared_distances=zeros(ntrain,1);
        %   for j=1:ntrain
        %       elementwise_delta=samples_train(j,:)-samples_test(j,:);
        %       squared_elementwise_delta=elementwise_delta.^2;
        %       squared_distance=sum(squared_elementwise_delta);
        %       squared_distances(j)=squared_distance;
        %   end

        [unused, i]=min(squared_distances);
        predicted(k)=targets_train(i);
        % <@@<
    end


text:
% AUTHOR:         Aaron Nicolson
% AFFILIATION:    Signal Processing Laboratory, Griffith University
%
% This Source Code Form is subject to the terms of the Mozilla Public
% License, v. 2.0. If a copy of the MPL was not distributed with this
% file, You can obtain one at http://mozilla.org/MPL/2.0/.

clear all; close all; clc;

%% GET MATLAB_FEAT REPOSITORY
addpath('~/Dropbox/GitHub/matlab_feat/feat')
addpath('./deepxi')

%% PARAMETERS
T_d = 32; % window duration (ms).
T_s = 16; % window shift (ms).
f_s = 16000; % sampling frequency (Hz).
s.N_d = round(f_s*T_d*0.001); % window duration (samples).
s.N_s = round(f_s*T_s*0.001); % window shift (samples).
s.f_s = f_s; % sampling frequency (Hz).
s.NFFT = 2^nextpow2(s.N_d); % frequency bins (samples).
d = s; x = s;
SNR_avg = -5:5:15; % SNR levels used to compute average SD level.

%% DIRECTORIES
s.dir = '/home/aaron/set/deep_xi_test_set/test_clean_speech';
d.dir = '/home/aaron/set/deep_xi_test_set/test_noise';
xi_hat_dir = input('xi_hat path:', 's');
xi_hat_dir_split = strsplit(xi_hat_dir, '/');
ver = [xi_hat_dir_split{end-2}, '_', xi_hat_dir_split{end-1}];

%% FILE LISTS
xi_hat_paths = dir([xi_hat_dir, '/*.mat']); % noise file paths.

results = MapNested();
noise_src_set = {};
SNR_set = {};
for i = 1:length(xi_hat_paths)

    load([xi_hat_paths(i).folder, '/', xi_hat_paths(i).name])

    if any(isnan(xi_hat(:))) || any(isinf(xi_hat(:)))
        error('NaN or Inf value in xi_hat.')
    end

    split_basename = strsplit(xi_hat_paths(i).name,'_');
    noise_src = split_basename{end-1};
    SNR = split_basename{end};
    clean_speech = extractBefore(xi_hat_paths(i).name, ['_', noise_src, '_', SNR]);
    SNR = SNR(1:end-6);

    s.wav = audioread([s.dir, '/', clean_speech, '_', noise_src, '.wav']); % clean speech.
    d.src = audioread([d.dir, '/', clean_speech, '_', noise_src, '.wav']); % noise.
    [~, d.wav] = add_noise(s.wav, d.src(1:length(s.wav)), str2double(SNR)); % noisy speech.

    s = analysis_stft(s, 'polar'); % clean-speech STMS.
    d = analysis_stft(d, 'polar'); % noise STMS.

    xi = (s.STMS.^2)./(d.STMS.^2); % instantaneous a priori SNR.

    xi_hat = xi_hat(1:size(xi, 1), :);
    D = spectral_distortion(xi, xi_hat);

    if any(isnan(D(:))) || any(isinf(D(:)))
        error('NaN or Inf value in D.')
    end

    if ~any(strcmp(SNR_set, SNR))
        SNR_set{end+1} = SNR;
    end

    if ~any(strcmp(noise_src_set, noise_src))
        noise_src_set{end+1} = noise_src;
    end

    if results.isKey(noise_src, SNR)
        results(noise_src, SNR) = [D; results(noise_src, SNR)];
    else
        results(noise_src, SNR) = D;
    end
    clc;
    fprintf('%.2f%%\n', 100*i/length(xi_hat_paths));
end

% there has to be a better way to do this.
for i=1:length(SNR_set)
    SNR_set{i} = str2num(SNR_set{i});
end
tmp = sort(cell2mat(SNR_set));
for i=1:length(SNR_set)
    SNR_set{i} = num2str(tmp(i));
end

res_dir = 'log/results/spectral_distortion_xi';
if ~exist(res_dir, 'dir')
    mkdir(res_dir)
end

fileID = fopen([res_dir, '/', ver, '.csv'],'w');
fprintf(fileID, 'noise, snr_db, SD\n');
avg = [];
for i=1:length(noise_src_set)
    for j=1:length(SNR_set)
        D = mean(results(noise_src_set{i}, SNR_set{j}));
        fprintf('%s, %s: %.2f.\n', noise_src_set{i}, SNR_set{j}, D);
        fprintf(fileID, '%s, %s, %.2f\n', noise_src_set{i}, SNR_set{j}, D);
        if ismember(j, SNR_avg)
            avg = [D; results(noise_src_set{i}, SNR_set{j})];
        end
    end
end
fclose(fileID);

avg_path = [res_dir, '/average.csv'];

if ~exist(avg_path, 'file')
    fileID = fopen(avg_path, 'w');
    fprintf(fileID, 'ver, SD\n');
    fclose(fileID);
end

fileID = fopen(avg_path, 'a');
fprintf(fileID, '%s, %.2f\n', ver, mean(avg));
fclose(fileID);

author: anicolson

repo: DeepXi

save_path: github-repos/MATLAB/anicolson-DeepXi


text:
% vgg_H_from_P_plane  Given 3D plane and camera, returns 4-by-3 matrix mapping image points to points on the plane.
%
% H = vgg_H_from_P_plane(A,P), where
%   A ... size (4,1), scene plane
%   P ... size (3,4), camera matrix
%   H ... size (4,3), matrix such that X = x*H, where x (size (3,1)) is an image point
%     and X is a scene point lying on the scene plane A.
%     In other words, P*H = s*eye(3), where s is a scalar factor.
%
% The overall sign of H is chosen such that s*A*vgg_wedge(P) > 0.

% T. Werner, Oct 2001

function H = vgg_H_from_P_plane(A,P)

A = A';
p1 = P(1,:);
p2 = P(2,:);
p3 = P(3,:);

H = [vgg_contreps(p2'*p3 - p3'*p2)*A ...
     vgg_contreps(p3'*p1 - p1'*p3)*A ...
     vgg_contreps(p1'*p2 - p2'*p1)*A];

return

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
function m_ungrid(goptn)
% M_UNGRID Removes a grid;
%          M_UNGRID deletes a map grid, but leaves any plotted
%          data.
%
%          M_UNGRID XXX
%          or
%          M_UNGRID('XXX')
%
%          can be used to remove other parameters plotted using an
%          M_XXX command (e.g. M_UNGRID COAST).

% Rich Pawlowicz (rich@ocgy.ubc.ca) 4/Apr/97
%
% 14/11/98 - Added possible option to remove other tagged items.
%
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.

if nargin==0
  mstr='m_grid_';
else
  mstr=['m_' lower(goptn)];
end

hh=get(gca,'children');

things=get(hh,'tag');

for i=1:length(hh)
  if ~isempty(things{i})  &&  strmatch(mstr,things{i})
   delete(hh(i));
 end
end

if strmatch('m_grid',mstr)
  set(gca,'visible','on');
end

author: CelsoReyes


text:
function x = gaussRnd(mu, Sigma, n)
% Generate samples from a Gaussian distribution.
% Input:
%   mu: d x 1 mean vector
%   Sigma: d x d covariance matrix
%   n: number of samples
% Outpet:
%   x: d x n generated sample x~Gauss(mu,Sigma)
% Written by Mo Chen (sth4nth@gmail.com).
if nargin == 2
    n = 1;
end
[V,err] = chol(Sigma);
if err ~= 0
    error('ERROR: sigma must be a symmetric positive definite matrix.');
end
x = V'*randn(size(V,1),n)+repmat(mu,1,n);

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function adj = triangulation_order3_adj_set ( node_num, ...
  tri_num, triangle_node, triangle_neighbor, adj_num, adj_col )

%*****************************************************************************80
%
%% TRIANGULATION_ORDER3_ADJ_SET sets adjacencies in a triangulation.
%
%  Discussion:
%
%    This routine is called to count the adjacencies, so that the
%    appropriate amount of memory can be set aside for storage when
%    the adjacency structure is created.
%
%    The triangulation is assumed to involve 3-node triangles.
%
%    Two nodes are "adjacent" if they are both nodes in some triangle.
%    Also, a node is considered to be adjacent to itself.
%
%    This routine can be used to create the compressed column storage
%    for a linear triangle finite element discretization of
%    Poisson's equation in two dimensions.
%
%  Diagram:
%
%       3
%    s  |\
%    i  | \
%    d  |  \
%    e  |   \  side 2
%       |    \
%    3  |     \
%       |      \
%       1-------2
%
%         side 1
%
%    The local node numbering
%
%
%   21-22-23-24-25
%    |\ |\ |\ |\ |
%    | \| \| \| \|
%   16-17-18-19-20
%    |\ |\ |\ |\ |
%    | \| \| \| \|
%   11-12-13-14-15
%    |\ |\ |\ |\ |
%    | \| \| \| \|
%    6--7--8--9-10
%    |\ |\ |\ |\ |
%    | \| \| \| \|
%    1--2--3--4--5
%
%    A sample grid
%
%
%    Below, we have a chart that summarizes the adjacency relationships
%    in the sample grid.  On the left, we list the node, and its neighbors,
%    with an asterisk to indicate the adjacency of the node to itself
%    (in some cases, you want to count this self adjacency and in some
%    you don't).  On the right, we list the number of adjancencies to
%    lower-indexed nodes, to the node itself, to higher-indexed nodes,
%    the total number of adjacencies for this node, and the location
%    of the first and last entries required to list this set of adjacencies
%    in a single list of all the adjacencies.
%
%    N   Adjacencies                Below  Self    Above  Total First  Last
%
%   --  -- -- -- -- -- -- --           --    --      --      --   ---     0
%    1:  *  2  6                        0     1       2       3     1     3
%    2:  1  *  3  6  7                  1     1       3       5     4     8
%    3:  2  *  4  7  8                  1     1       3       5     9    13
%    4:  3  *  5  8  9                  1     1       3       5    14    18
%    5:  4  *  9 10                     1     1       2       4    19    22
%    6:  1  2  *  7 11                  2     1       2       5    23    27
%    7:  2  3  6  *  8 11 12            3     1       3       7    28    34
%    8:  3  4  7  *  9 12 13            3     1       3       7    35    41
%    9:  4  5  8  * 10 13 14            3     1       3       7    42    48
%   10:  5  9  * 14 15                  2     1       2       5    49    53
%   11:  6  7  * 12 16                  2     1       2       5    54    58
%   12:  7  8 11  * 13 16 17            3     1       3       7    59    65
%   13:  8  9 12  * 14 17 18            3     1       3       7    66    72
%   14:  9 10 13  * 15 18 19            3     1       3       7    73    79
%   15: 10 14  * 19 20                  2     1       2       5    80    84
%   16: 11 12  * 17 21                  2     1       2       5    85    89
%   17: 12 13 16  * 18 21 22            3     1       3       7    90    96
%   18: 13 14 17  * 19 22 23            3     1       3       7    97   103
%   19: 14 15 18  * 20 23 24            3     1       3       7   104   110
%   20: 15 19  * 24 25                  2     1       2       5   111   115
%   21: 16 17  * 22                     2     1       1       4   116   119
%   22: 17 18 21  * 23                  3     1       1       5   120   124
%   23: 18 19 22  * 24                  3     1       1       5   125   129
%   24: 19 20 23  * 25                  3     1       1       5   130   134
%   25: 20 24  *                        2     1       0       3   135   137
%   --  -- -- -- -- -- -- --           --    --      --      --   138   ---
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer TRI_NUM, the number of triangles.
%
%    Input, integer TRIANGLE_NODE(3,TRI_NUM), lists the nodes that
%    make up each triangle in counterclockwise order.
%
%    Input, integer TRIANGLE_NEIGHBOR(3,TRI_NUM), for each side of
%    a triangle, lists the neighboring triangle, or -1 if there is
%    no neighbor.
%
%    Input, integer ADJ_NUM, the number of adjacencies.
%
%    Input, integer ADJ_COL(NODE_NUM+1).  Information about column J is stored
%    in entries ADJ_COL(J) through ADJ_COL(J+1)-1 of ADJ.
%
%    Output, integer ADJ(ADJ_NUM), the adjacency information.
%
  triangle_order = 3;
  adj(1:adj_num) = -1;
  adj_copy(1:node_num) = adj_col(1:node_num);
%
%  Set every node to be adjacent to itself.
%
  for node = 1 : node_num
    adj(adj_copy(node)) = node;
    adj_copy(node) = adj_copy(node) + 1;
  end
%
%  Examine each triangle.
%
  for triangle = 1 : tri_num

    n1 = triangle_node(1,triangle);
    n2 = triangle_node(2,triangle);
    n3 = triangle_node(3,triangle);
%
%  Add edge (1,2) if this is the first occurrence,
%  that is, if the edge (1,2) is on a boundary (TRIANGLE2 <= 0)
%  or if this triangle is the first of the pair in which the edge
%  occurs (TRIANGLE < TRIANGLE2).
%
    triangle2 = triangle_neighbor(1,triangle);

    if ( triangle2 < 0 | triangle < triangle2 )
      adj(adj_copy(n1)) = n2;
      adj_copy(n1) = adj_copy(n1) + 1;
      adj(adj_copy(n2)) = n1;
      adj_copy(n2) = adj_copy(n2) + 1;
    end
%
%  Add edge (2,3).
%
    triangle2 = triangle_neighbor(2,triangle);

    if ( triangle2 < 0 | triangle < triangle2 )
      adj(adj_copy(n2)) = n3;
      adj_copy(n2) = adj_copy(n2) + 1;
      adj(adj_copy(n3)) = n2;
      adj_copy(n3) = adj_copy(n3) + 1;
    end
%
%  Add edge (3,1).
%
    triangle2 = triangle_neighbor(3,triangle);

    if ( triangle2 < 0 | triangle < triangle2 )
      adj(adj_copy(n1)) = n3;
      adj_copy(n1) = adj_copy(n1) + 1;
      adj(adj_copy(n3)) = n1;
      adj_copy(n3) = adj_copy(n3) + 1;
    end

  end
%
%  Ascending sort the entries for each node.
%
  for node = 1 : node_num
    k1 = adj_col(node);
    k2 = adj_col(node+1)-1;
    adj(k1:k2) = i4vec_sort_heap_a ( k2+1-k1, adj(k1:k2) );
  end

  return
end


text:
function output = callclp(interfacedata)

% Retrieve needed data
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
ub      = interfacedata.ub;
lb      = interfacedata.lb;
Q       = interfacedata.Q;

n = length(c);

if options.showprogress;showprogress(['Calling CLP'],options.showprogress);end

if K.f>0
    Aeq = -F_struc(1:K.f,2:end);
    beq = F_struc(1:K.f,1);
else
    Aeq = [];
    beq = [];
end

if any(K.l)
    A = -F_struc(1+K.f:end,2:end);
    b = F_struc(1+K.f:end,1);
else
    A = [];
    b = [];
end

% Fix for bug in mexclp or clp
fixed = 0;
if size(A,1)==0 & size(Aeq,2)==0
    A = [ones(1,length(c))];
    if ~isempty(ub)
        A(isinf(ub)) = 0;
        dummy = ub;
        dummy(isinf(ub)) = 0;
        b = sum(dummy);
    else
        b = 1e8;
    end
    fixed = 1;
    options.saveduals = 0;
end
% lb(lb==-inf)=-1e12
% ub(ub==inf)=1e12

ops = options.clp;
ops.verbose = options.verbose;

if options.savedebug
    save clpdebug
end

% Call mex-interfacec
solvertime = tic;
[x,lambda,problem] = clp(2*Q,c,A,b,Aeq,beq,lb,ub,ops);%,interfacedata.integer_variables);
solvertime = toc(solvertime);

if options.saveduals
    D_struc = -lambda;
else
    D_struc = [];
end

% Save all data sent to solver?
if options.savesolverinput
	solverinput = [];
else
	solverinput = [];
end

% Save all data from the solver?
if options.savesolveroutput
	solveroutput = [];
else
	solveroutput = [];
end

% Standard interface
output = createOutputStructure(x,D_struc,[],problem,interfacedata.solver.tag,solverinput,solveroutput,solvertime);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function X=sos(X,r)
%SOS Declare sum-of-squares structure
%
% F = sos(p,r)
%
% Input
%  p : SDPVAR object
%  r : Desired rank (optional)
% Output
%  F : Constraint object
%
% Example:
%  Typical usage is
%
%   F = sos(p)
%
%  An experimental feature is to search for
%  low rank decompositions. To search for a
%  decomposition using at most 3 terms, use
%  a second argument
%
%   F = sos(p,3)
%
%  Note that his feature requires the solver LMIRANK.
%
% See also SOSEX

if ~issymmetric(X)% X.n>1 | X.m>1
    error('SOS can only be applied to symmetric polynomial matrices');
end
if nargin<2
    r = inf;
end
X.typeflag = 11;
X.extra.sosid = yalmip('sosid');
X.extra.rank = r;

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [fx,dfdx,dfdp] = f_HRF2(Xt,P,ut,in)
% Balloon (HRF) model evolution function in log-space
% function [fx,dfdx,dfdp] = f_HRF2(Xt,P,ut,in)
% This function evaluates the evolution function derived from the balloon
% model for the hemodynamic response function. It can be called in two
% ways: (i) as a "stand-alone" evolution function, whereby the system's
% states are hemodynamic states of the balloon model, or (ii) as a
% generalized observation function, where the real system's states are the
% neuronal states of a DCM for fMRI model. Note that the hemodynamic states
% are in log-space, for positivity constrints.

deltat = in.deltat;
n = size(Xt,1);

% Get parameters
[E0,V0,tau0,kaf,kas,epsilon,alpha] = BOLD_parameters;
if isfield(in,'fullDCM') && in.fullDCM
    nreg = n./5;
    ind1 = in.ind1;
    ind2 = in.ind2;
    ind3 = in.ind3;
    ind4 = in.ind4;
    n1 = in.n1;
    n2 = in.n2;
    n3 = in.n3;
    n4 = in.n4;
    try, n5=in.n5; catch, n5=[];end
    epsilon = 1;
    alpha = alpha.*exp(P(in.ind5));
else
    nreg = n./4;
    ind1 = 1;
    ind2 = 2;
    ind3 = 3;
    ind4 = 4;
    n1 = 1;
    n2 = 2;
    n3 = 3;
    n4 = 4;
    epsilon = epsilon.*exp(P(5));
    alpha = alpha.*exp(P(6));
end
% [E0,dsdp] = sigm(P(ind1)-0.6633,struct('beta',1,'G0',1,'INV',0));
% E0 = E0(:);
E0 = 1./(1+exp(-(P(ind1)-0.6633)));
dsdp = E0.*(1-E0);
tau0 = tau0.*exp(P(ind2));
kaf = kaf.*exp(P(ind3));
kas = kas.*exp(P(ind4));

if isfield(in,'xshift') % for numerical stability
    xshift = in.xshift;
else
    xshift = 0;
end

% hemodynamic states ...
if isfield(in,'linearized') && in.linearized
    x1 = zeros(nreg,1);
    x2 = ones(nreg,1);
    x3 = ones(nreg,1);
    x4 = ones(nreg,1);
else
    % vasodilatory signal s(t)
    x1 = Xt(n1,:);
    % blood inflow f(t)
    if isfield(in,'logx2') && ~in.logx2
        x2 = Xt(n2,:) + 1; % deviation to steady-state!
    else
        x2 = exp(Xt(n2,:)) + xshift;
    end
    % blood volume v(t)
    x3 = exp(Xt(n3,:)) + 0.*xshift;
    % dHb content q(t)
    x4 = exp(Xt(n4,:)) + 0.*xshift;
end
% blood outflow
fv = x3.^(1./alpha);
% d[blood flow]/dXt(3)
dfvdx = (1./alpha).*fv;
% oxygen extraction
ff = (1-(1-E0).^(1./x2))./E0;
% d[O2 extraction]/dXt(2)
if isfield(in,'logx2') && ~in.logx2
    dffdx = log(1-E0).*(1-E0).^(1./x2)./(E0.*x2.^2);
else
    dffdx = log(1-E0).*(1-E0).^(1./x2)./(E0.*x2);
end
% ... and flow field, derivatives, etc...
f = zeros(n,1);
J = zeros(n,n);
dfdp = zeros(size(P,1),n);

% Evaluate flow field
f(n1) = (epsilon.*ut - kas.*x1 - kaf.*(x2 - 1));
if isfield(in,'logx2') && ~in.logx2
    f(n2) = x1;
else
    f(n2) = x1./x2;
end
f(n3) = (x2 - fv)./(tau0.*x3);
f(n4) = (x2.*ff./x4 - fv./x3)./tau0;

% Evaluate jacobian and gradients wrt parameters
for i=1:nreg

    if isfield(in,'logx2') && ~in.logx2
        J(n1(i),n1(i):n1(i)+3) = [ -kas(i) , 1, 0, 0 ];
        J(n2(i),n1(i):n1(i)+3) = [ -kaf(i), 0, ...
            1./(tau0(i).*x3(i)), ...
            (ff(i)+x2(i).*dffdx(i))./(x4(i).*tau0(i))];
    else
        J(n1(i),n1(i):n1(i)+3) = [ -kas(i) , 1./x2(i), 0, 0 ];
        J(n2(i),n1(i):n1(i)+3) = [ -kaf(i).*x2(i), -x1(i)./x2(i), ...
            x2(i)./(tau0(i).*x3(i)), ...
            x2(i).*(ff(i)+dffdx(i))./(tau0(i).*x4(i))];
    end
    J(n3(i),n1(i):n1(i)+3) = [ 0, 0,...
        -f(n3(i)) - dfvdx(i)./(tau0(i).*x3(i)) , ...
        (fv(i)-dfvdx(i))./(tau0(i).*x3(i))];
    J(n4(i),n1(i):n1(i)+3) = [ 0, 0, 0, ...
        -(x2(i).*ff(i))./(tau0(i).*x4(i)) ];

    if isfield(in,'linearized') && in.linearized

        tmp = log(1-E0(i))./E0(i);

        dfdp(ind1(i),n1(i):n1(i)+3) = [0,0,0,...
            -dsdp(i).*(1+tmp).*Xt(n2(i))./(tau0(i).*E0(i))];
        dfdp(ind2(i),n1(i):n1(i)+3) = [0,0,...
            (-Xt(n2(i))+Xt(n3(i))./alpha(i))./tau0(i),...
            -((ff(i)+dffdx(i)).*Xt(n2(i)) + (fv(i)-dfvdx(i)).*Xt(n3(i)) ...
            - ff(i).*Xt(n4(i)))./tau0(i)];
        dfdp(ind3(i),n1(i):n1(i)+3) = kaf(i).*[-Xt(n2(i)),0,0,0];
        dfdp(ind4(i),n1(i):n1(i)+3) = kas(i).*[-Xt(n1(i)),0,0,0];

        if isfield(in,'fullDCM') && in.fullDCM
            if  ~isempty(n5)
                J(n5(i),n1(i)) = epsilon;
            end
            dfdp(in.ind5(i),n1(i):n1(i)+3) = [0,0,...
                Xt(n3(i))./(tau0(i).*alpha(i)),...
                Xt(n3(i))./(tau0(i).*alpha(i))];
        else
            dfdp(5,:) = epsilon.*[ut,0,0,0];
            dfdp(6,:) = [0,0,...
                Xt(n3(i))./(tau0(i).*alpha(i)),...
                Xt(n3(i))./(tau0(i).*alpha(i))];
        end

    else

        % gradient wrt parameters
        dfdp(ind1(i),n1(i):n1(i)+3) = [0,0,0,...
            (((1-E0(i)).^(-1+1./x2(i)))-x2(i).*ff(i))...
            .*dsdp(i)./(tau0(i).*x4(i).*E0(i))];
        dfdp(ind2(i),n1(i):n1(i)+3) = [0,0,...
            -(x2(i) - fv(i))./(tau0(i).*x3(i)),...
            -(x2(i).*ff(i)./x4(i) - fv(i)./x3(i))./tau0(i)];
        dfdp(ind3(i),n1(i):n1(i)+3) = kaf(i).*[-x2(i)+1,0,0,0];
        dfdp(ind4(i),n1(i):n1(i)+3) = kas(i).*[-x1(i),0,0,0];

        % complement gradients if used for full DCM model
        if isfield(in,'fullDCM') && in.fullDCM
            if  ~isempty(n5)
                J(n5(i),n1(i)) = epsilon;
            end
            dfdp(in.ind5(i),n1(i):n1(i)+3) = [0,0,...
                log(x3(i)).*fv(i)./(tau0(i).*x3(i).*alpha(i)),...
                log(x3(i)).*fv(i)./(tau0(i).*x3(i).*alpha(i))];
        else
            dfdp(5,:) = epsilon.*[ut,0,0,0];
            dfdp(6,:) = [0,0,...
                log(x3).*fv./(tau0.*x3.*alpha),...
                log(x3).*fv./(tau0.*x3.*alpha)];
        end

    end

end

% Apply Euler discretization
if isfield(in,'linearized') && in.linearized
     if isfield(in,'fullDCM') && in.fullDCM && ~isempty(n5)
         Cu = 0;
     else
         Cu = zeros(n,1);
         Cu(n1) = epsilon.*ut;
     end
    fx = Xt + deltat.*(J'*Xt + Cu);
else
    fx = Xt + deltat.*f;
end
dfdp = deltat.*dfdp;
dfdx = eye(n) + deltat.*J;

author: MBB-team

repo: VBA-toolbox


text:
function id = defaultId(stack)

% DEFAULTID returns a string that can serve as warning or error identifier,
% for example 'FieldTip:ft_read_header:line345'.
%
% See also WARNING, ERROR, FT_NOTICE, FT_INFO, FT_DEBUG

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

if nargin==1
  % stack is already in input, and assumed to be properly pruned
else
  stack = dbstack('-completenames');

  % stack(1).name is defaultId
  % stack(2).name is ft_notification
  % stack(3).name is ft_error, ft_warning, ft_notice, ft_info or ft_debug

  % remove the functions that pertain to the notification system itself
  stack = stack(4:end);

  % remove the non-FieldTrip functions and scripts, these should not be part of the message identifier
  [v, p] = ft_version;
  keep   = startsWith({stack.file}, p);
  stack  = stack(keep);
end

switch numel(stack)
  case 0
    % it is called from the command line
    id = 'FieldTrip:commandline';
  case 1
    id = sprintf('FieldTrip:%s:line%d', stack(1).name, stack(1).line);
  case 2
    id = sprintf('FieldTrip:%s:%s:line%d', stack(2).name, stack(1).name, stack(1).line);
  otherwise
    % it is called from within a function
    name = sprintf('%s:', stack(end:-1:1).name); % this creates something like fun1:fun2:fun3:
    id   = sprintf('FieldTrip:%sline%d', name, stack(1).line);
end

% slashes occur when using nested functions, but are not allowed in the identifier
id(id=='/') = ':';
id(id=='\') = ':';

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% example_nufft1_antenna.m
%
% This m-file is an example of applying the NUFFT method to compute
% the far-field beam pattern of a 2D array of nonuniformly spaced
% (point?) antenna elements.
% In the nomenclature of the 1999 SIAM J Sci. Comput. paper by Nguyen and Liu,
% this is "Problem 1."
%
% This provides an example of how to compute "the Fourier transform"
% of nonuniformly spaced data.
%
% The user should think very carefully about whether it is truly meaningful
% to compute "the FT" of unequally spaced data first.  It certainly is
% meaningful for attenna beam patterns, but may not always be what is really
% useful for all applications with nonuniformly sampled data.
%
% Copyright 2007, Jeff Fessler, University of Michigan

%
% antenna element locations in 2d plane
%
clf, pl = @(p) subplot(220+p);
for choice=1:2
	tmp = [0:40]'/41 * 2 * pi;
	yc = sin(choice*tmp); % funny bowtie pattern for illustration
	xc = cos(tmp); clear tmp
	pl(choice), plot(xc, yc, 'o'), axis square
	xlabel 'x', ylabel 'x', title 'antenna locations'

	% create NUFFT structure
	N = [1 1]*2^8;
	J = [5 5];	% interpolation neighborhood
	K = N*2;	% two-times oversampling
	om = [xc yc];	% 'frequencies' are locations here!

	% the following line probably needs to be changed
	% to get the proper scaling/dimensions in the pattern
	% but for now i just make it fill up [-pi/20,pi/20]
	% in hopes of getting a nice 'picture' of pattern
	om = (pi/20) * om / max(om(:));
	st = nufft_init(om, N, J, K, N/2, 'minmax:kb');

	weights = ones(size(xc)); % equal weights on each element; could change

	% call the *adjoint* NUFFT; this is what does "the FT of unequal data"
	pattern = nufft_adj(weights, st);

	pl(2+choice)
	im(pattern, 'pattern')
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function output = linear_transform_animation(data,A)
% data = N x 2 matrix
% A = 2 x 2 matrix
% output = N x 2 matrix

dot_colors = jet(size(data,1));

output = transpose(A*data');

minmin = min(min(data));
maxmax = max(max(data));

xy_min = -max(abs(minmin),abs(maxmax));
xy_max = max(abs(minmin),abs(maxmax));

figure;
set(gcf,'position',[316.2000  215.4000  904.8000  420.0000]);
subplot(1,2,1);
scatter(data(:,1),data(:,2),30,dot_colors,'filled')
grid on;
xlim([xy_min, xy_max]); ylim([xy_min, xy_max]);
set(gcf,'color','white');

subplot(1,2,2);
scatter(data(:,1),data(:,2),30,dot_colors,'filled')
grid on;
xlim([xy_min, xy_max]); ylim([xy_min, xy_max]);
set(gcf,'color','white');

n_steps = 100;
for i_steps = 1:n_steps
    step_mtx = (A-eye(2))/n_steps*i_steps;
    new_xy = (eye(2)+step_mtx)*[data(:,1), data(:,2)]';
    scatter(new_xy(1,:), new_xy(2,:),30,dot_colors,'filled')
    grid on;
    xlim([xy_min, xy_max]); ylim([xy_min, xy_max]);
    set(gcf,'color','white')
    pause(0.001);
end

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
function [ opts ] = PrepareData_Char_RNN( opts )

max_char=67;

x = textread('train_x.txt');
x(:,end)=[];
opts.train = zeros(max_char, size(x,1),size(x,2));
Index=x(:)'+1+max_char*[(0:numel(x)-1)];
opts.train(Index)=1;

x = textread('train_y.txt');
opts.train_labels=x+1;

opts.n_train=size(opts.train_labels,1);

x = textread('test_x.txt');
x(:,end)=[];
opts.test = zeros(max_char, size(x ,1),size(x ,2));
Index=x (:)'+1+max_char*[(0:numel(x )-1)];
opts.test(Index)=1;

x = textread('test_y.txt');
opts.test_labels=x+1;

opts.n_test=size(opts.test_labels,1);

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
function determ = daub12_determinant ( n )

%*****************************************************************************80
%
%% DAUB12_DETERMINANT returns the determinant of the DAUB12 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Output, real DETERM, the determinant.
%
  determ = - 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
 function y = ir_dct2(x)
%function y = ir_dct2(x)
%|
%| 2D DCT via separable 1D processing for many 2D images
%| because matlab dct2() only handles a single 2D image.
%| in
%| x [nx ny npatch]
%| out
%| y [nx ny npatch]
%|
%| 2017 Jeff Fessler, University of Michigan

if nargin ~= 1, ir_usage, end

[nx, ny, np] = size(x);
y = permute(reshape(...
	dct(reshape(permute(reshape(...
	dct(reshape(x, nx, [])), nx, ny, np), ...
	[2 1 3]), ny, [])), ny, nx, np), [2 1 3]);


text:
function sin_degree_values_test ( )

%*****************************************************************************80
%
%% SIN_DEGREE_VALUES_TEST demonstrates the use of SIN_DEGREE_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 January 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SIN_DEGREE_VALUES_TEST:\n' );
  fprintf ( 1, '  SIN_DEGREE_VALUES stores values of \n' );
  fprintf ( 1, '  the sine function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X           FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = sin_degree_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end