text:
function [] = test_l1_logistic_regression()

    clc;
    clear;
    close all;

    %rng('default');

    %% Set algorithms
    if 0
        algorithms = sgd_solver_list('ALL');
    else
        algorithms = {'SVRG', 'SAG', 'Adam', 'Reg-oBFGS-Lim','APG-BKT', 'APG-TFOCS-BKT'};
    end

    %% prepare dataset
    if 1
        % generate synthtic data
        d = 100;
        n = 1000;
        data = logistic_regression_data_generator(n, d);
        x_train = data.x_train;
        y_train = data.y_train;
        x_test = data.x_test;
        y_test = data.y_test;
        d = size(x_train,1);
        %w_opt = data.w_opt;
        lambda = 0.1;
    else
        % load pre-created synthetic data
        data = importdata('../data/logistic_regression/data_100d_10000.mat');
        x_train = data.x_train;
        y_train = data.y_train;
        x_test = data.x_test;
        y_test = data.y_test;
        d = size(x_train,1);
        n = length(y_train);
        w_opt = data.w_star;
        lambda = data.lambda;
    end

    %% define problem definitions
    problem = l1_logistic_regression(x_train, y_train, x_test, y_test, lambda);

    %% initialize
    w_init = rand(d,1);
    w_list = cell(length(algorithms),1);
    info_list = cell(length(algorithms),1);

    %% calculate solution
    %if norm(w_opt)
    %else
        % calculate solution
        options.w_init = w_init;
        options.max_iter = 100;
        w_opt = problem.calc_solution(options, 'ag');
    %end
    f_opt = problem.cost(w_opt);
    fprintf('f_opt: %.24e\n', f_opt);

    %% calculate solution
    f_opt = problem.cost(w_opt);
    fprintf('f_opt: %.24e\n', f_opt);

    %% perform algorithms
    for alg_idx=1:length(algorithms)
        fprintf('\n\n### [%02d] %s ###\n\n', alg_idx, algorithms{alg_idx});

        clear options;
        % general options for optimization algorithms
        options.w_init = w_init;
        options.tol_gnorm = 1e-10;
        options.max_iter = 100;
        options.max_epoch = 100;
        options.verbose = true;
        options.f_opt = f_opt;
        batch_size = 10;

        switch algorithms{alg_idx}
            case {'PG-BKT'}

                options.step_alg = 'backtracking';
                options.step_init_alg = 'bb_init';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'PG-TFOCS-BKT'}

                options.step_alg = 'tfocs_backtracking';
                options.step_init_alg = 'bb_init';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'APG-BKT'}

                options.step_alg = 'backtracking';
                options.step_init_alg = 'bb_init';
                [w_list{alg_idx}, info_list{alg_idx}] = ag(problem, options);

            case {'APG-TFOCS-BKT'}

                options.step_alg = 'tfocs_backtracking';
                options.step_init_alg = 'bb_init';
                [w_list{alg_idx}, info_list{alg_idx}] = ag(problem, options);

            case {'L-BFGS-BKT'}

                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            case {'L-BFGS-WOLFE'}

                options.step_alg = 'strong_wolfe';
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            case {'L-BFGS-TFOCS'}

                options.step_alg = 'tfocs_backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            case {'BFGS-TFOCS'}

                options.step_alg = 'tfocs_backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'Newton-CHOLESKY'}

                options.sub_mode = 'CHOLESKY';
                options.step_alg = 'backtracking';
                %options.step_alg = 'tfocs_backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = newton(problem, options);

            case {'NCG-BKT'}

                options.sub_mode = 'STANDARD';
                options.step_alg = 'backtracking';
                %options.step_alg = 'tfocs_backtracking';
                %options.beta_alg = 'PR';
                [w_list{alg_idx}, info_list{alg_idx}] = ncg(problem, options);

            case {'SGD'}

                options.batch_size = batch_size;
                options.step = 0.001 * options.batch_size;
                %options.step_alg = 'decay';
                options.step_alg = 'fix';

                [w_list{alg_idx}, info_list{alg_idx}] = sgd(problem, options);

            % Variance reduction (VR) varitns
            case {'SVRG'}

                options.batch_size = batch_size;
                options.step_init = 0.0001 * options.batch_size;
                options.step_alg = 'fix';

                [w_list{alg_idx}, info_list{alg_idx}] = svrg(problem, options);

            case {'SAG'}

                options.batch_size = batch_size;
                %options.step_init = 0.00005 * options.batch_size;
                options.step_init = 0.0001 * options.batch_size;
                options.step_alg = 'fix';
                options.sub_mode = 'SAG';

                [w_list{alg_idx}, info_list{alg_idx}] = sag(problem, options);

            case {'SAGA'}

                options.batch_size = batch_size;
                %options.step_init = 0.00005 * options.batch_size;
                options.step_init = 0.000001 * options.batch_size;
                options.step_alg = 'fix';
                options.sub_mode = 'SAGA';

                [w_list{alg_idx}, info_list{alg_idx}] = sag(problem, options);

            % AdaGrad variants
            case {'AdaGrad'}

                options.batch_size = batch_size;
                options.step_init = 0.02 * options.batch_size;
                options.step_alg = 'fix';
                options.epsilon = 0.00001;
                options.sub_mode = 'AdaGrad';

                [w_list{alg_idx}, info_list{alg_idx}] = adagrad(problem, options);

            case {'RMSProp'}

                options.batch_size = batch_size;
                options.step_init = 0.00001 * options.batch_size;
                options.step_alg = 'fix';
                options.epsilon = 0.00001;
                options.sub_mode = 'RMSProp';
                options.beta = 0.9;

                [w_list{alg_idx}, info_list{alg_idx}] = adagrad(problem, options);

            case {'AdaDelta'}

                options.batch_size = batch_size;
                options.step_init = 0.01 * options.batch_size;
                options.step_alg = 'fix';
                options.epsilon = 0.00001;

                options.sub_mode = 'AdaDelta';
                options.beta = 0.9;

                [w_list{alg_idx}, info_list{alg_idx}] = adagrad(problem, options);

            case {'Adam'}

                options.batch_size = batch_size;
                options.step_init = 0.00001 * options.batch_size;
                options.step_alg = 'fix';
                options.sub_mode = 'Adam';
                options.beta1 = 0.8;
                options.beta2 = 0.999;
                options.epsilon = 0.00001;

                [w_list{alg_idx}, info_list{alg_idx}] = adam(problem, options);

            case {'AdaMax'}

                options.batch_size = batch_size;
                options.step_init = 0.00001 * options.batch_size;
                options.step_alg = 'fix';
                options.sub_mode = 'AdaMax';
                options.beta1 = 0.8;
                options.beta2 = 0.999;
                options.epsilon = 0.00001;

                [w_list{alg_idx}, info_list{alg_idx}] = adam(problem, options);

            % Stochastic Quasi-Newton variants
            case {'SQN'}

                options.batch_size = batch_size;
                options.batch_hess_size = batch_size * 20;
                options.step_init = 0.001 * options.batch_size;
                options.step_alg = 'fix';
                options.sub_mode = 'SQN';
                options.L = 20;
                options.mem_size = 20;

                [w_list{alg_idx}, info_list{alg_idx}] = slbfgs(problem, options);

            case {'SVRG-SQN'}

                options.batch_size = batch_size;
                options.batch_hess_size = batch_size * 20;
                options.step_init = 0.01 * options.batch_size;
                options.step_alg = 'fix';
                options.sub_mode = 'SVRG-SQN';
                options.L = 20;
                options.mem_size = 20;

                [w_list{alg_idx}, info_list{alg_idx}] = slbfgs(problem, options);

            case {'SVRG-LBFGS'}

                options.batch_size = batch_size;
                options.batch_hess_size = batch_size * 20;
                options.step_init = 0.01 * options.batch_size;
                options.step_alg = 'fix';
                options.sub_mode = 'SVRG-LBFGS';
                options.mem_size = 20;

                [w_list{alg_idx}, info_list{alg_idx}] = slbfgs(problem, options);

            case {'SS-SVRG'}

                options.batch_size = batch_size;
                options.batch_hess_size = batch_size * 20;
                options.step_init = 0.0005 * options.batch_size;
                options.step_alg = 'fix';
                r = d-1;
                if r < 1
                    r = 1;
                end
                options.r = r;

                [w_list{alg_idx}, info_list{alg_idx}] = subsamp_svrg(problem, options);

            case {'oBFGS-Inf'}

                options.batch_size = batch_size;
                options.step_init = 0.0001 * options.batch_size;
                options.step_alg = 'fix';
                options.sub_mode = 'Inf-mem';
                options.regularized = false;

                [w_list{alg_idx}, info_list{alg_idx}] = obfgs(problem, options);

            case {'oBFGS-Lim'}

                options.batch_size = batch_size;
                options.step_init = 0.00001 * options.batch_size;
                options.step_alg = 'fix';
                options.sub_mode = 'Lim-mem';
                options.mem_size = 20;
                options.regularized = false;

                [w_list{alg_idx}, info_list{alg_idx}] = obfgs(problem, options);

            case {'Reg-oBFGS-Inf'}

                options.batch_size = batch_size;
                options.step_init = 0.0001 * options.batch_size;
                options.step_alg = 'fix';
                options.sub_mode = 'Inf-mem';
                options.regularized = true;
                options.delta = 0.1;

                [w_list{alg_idx}, info_list{alg_idx}] = obfgs(problem, options);

            case {'Reg-oBFGS-Lim'}

                options.batch_size = batch_size;
                options.step_init = 0.0001 * options.batch_size;
                options.step_alg = 'fix';
                options.sub_mode = 'Lim-mem';
                options.mem_size = 20;
                options.regularized = true;
                options.delta = 0.1;

                [w_list{alg_idx}, info_list{alg_idx}] = obfgs(problem, options);

            case {'Damp-oBFGS-Inf'} % SDBFGS

                options.batch_size = batch_size;
                options.step_init = 0.0001 * options.batch_size;
                options.step_alg = 'fix';
                options.sub_mode = 'Inf-mem';
                options.regularized = true;
                options.delta = 0.1;
                options.damped = true;

                [w_list{alg_idx}, info_list{alg_idx}] = obfgs(problem, options);

            case {'Damp-oBFGS-Lim'}

                options.batch_size = batch_size;
                options.step_init = 0.01 * options.batch_size;
                options.step_alg = 'fix';
                options.sub_mode = 'Lim-mem';
                options.regularized = true;
                options.delta = 0.1;
                options.damped = true;

                [w_list{alg_idx}, info_list{alg_idx}] = obfgs(problem, options);

           case {'IQN'}

                options.w_init = w_init;
                options.step_init = 1;
                options.step_alg = 'fix';

                [w_list{alg_idx}, info_list{alg_idx}] = iqn(problem, options);

            otherwise
                warn_str = [algorithms{alg_idx}, ' is not supported.'];
                warning(warn_str);
                w_list{alg_idx} = '';
                info_list{alg_idx} = '';
        end

    end

    fprintf('\n\n');

    %% plot all
    close all;

    % display iter vs cost/gnorm
    display_graph('iter','cost', algorithms, w_list, info_list);
    % display iter vs. l1 norm, i.e. the toral number of non-zero elements
    display_graph('iter','reg', algorithms, w_list, info_list);

    display_graph('grad_calc_count','optimality_gap', algorithms, w_list, info_list);

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function CM=ImaGIN_connectivity_matrix(N)
% -=============================================================================
% This function is part of the ImaGIN software:
% https://f-tract.eu/
%
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE AUTHORS
% DO NOT ASSUME ANY LIABILITY OR RESPONSIBILITY FOR ITS USE IN ANY CONTEXT.
%
% Copyright (c) 2000-2018 Inserm U1216
% =============================================================================-
%
% Authors: Olivier David

ncouple	= N*(N-1)/2;
if ncouple>0
    CM	= zeros(2,ncouple);
    cou	= 0;
    for ii	= 1:N-1
        j	= ii+1;
        cou	= cou+1;
        CM(1,cou)	= ii;
        CM(2,cou)	= j;
        while j < N
            j	= j+1;
            cou	= cou+1;
            CM(1,cou)	= ii;
            CM(2,cou)	= j;
        end
    end
else
    CM=[];
end
return

author: brainstorm-tools


text:
% Small technical example.
% Performs a check on the computation of certain central moments
% on a quincunx grid.
%
disp('Small technical example.');
disp('Performs a check on the computation of certain central moments');
disp('on a quincunx grid.');
disp('FOR MORE INFORMATION:  help Q0011mupq');
disp(' ');
disp('See also the report http://repository.cwi.nl:8888/cwi_repository/docs/IV/04/04178D.pdf');
disp('Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>');
disp(' (C) 1998-2006 Stichting CWI, Amsterdam, The Netherlands');
disp(' ');
%---INSERT YOUR IMAGE HERE-----------------------------------------------------
if exist('imread','file') == 2
  Orig = double(imread('zenithgray.TIF','tiff'));
else
  load zenithgray; Orig = zenithgray; clear zenithgray;
end
%
disp([' Dimensions ' int2str( size(Orig) )]);
%---EXTRACT A QUINCUNX GRIDFUNCTION-----
F00 = getcolor00(Orig);
F11 = getcolor11(Orig);
disp('The dimensions of the extracted quincunx gridfunction read as follows');
disp([int2str(size(F00)) ' and ' int2str(size(F11))]);
%
%---COMPUTE MASS------------------------
mu00 = Q0011mupq(F00, F11, 0, 0);
disp([' Mass  ' num2str( mu00, '%+12.4e')]);
%---NICE CHECK: OUTCOME SHOULD VANISH---
mu10 = Q0011mupq(F00, F11, 1, 0);
disp(' NICE CHECK: First order central moment (x-dir) should vanish');
disp([' reads ' num2str( mu10, '%+12.4e')]);
%---NICE CHECK: OUTCOME SHOULD VANISH---
mu01 = Q0011mupq(F00, F11, 0, 1);
disp(' NICE CHECK: First order central moment (y-dir) should vanish');
disp([' reads ' num2str( mu01, '%+12.4e')]);

author: Sable

repo: mcbench-benchmarks


text:
function varargout = par_testworker(address,func,varargin)
% Test a remote worker process by forwarding a task.
% Result = par_testworker(Address,Function,Arguments...)
%
% In:
%   Address : network address of the worker (Host:Port)
%
%   Function : handle of the function to evaluate
%
%   Arguments... : optional arguments to the function
%
% Out:
%   Result : output of the function
%
% Examples:
% 	par_testworker('192.168.0.1:23457',@sin,1000)
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2013-07-18

% Copyright (C) Christian Kothe, SCCN, 2013, christian@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify it under the terms of the GNU
% General Public License as published by the Free Software Foundation; either version 2 of the
% License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without
% even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License along with this program; if not,
% write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
% USA

import java.io.*
import java.net.*
import java.lang.*

% some parameters
taskid = round(rand*100000);
disp(['assigning task id ' num2str(taskid)]);
receive_buffer_size = 64000;
receive_backlog = 10;
timeout_connect_ms = 3000;
timeout_send_ms = 5000;
timeout_recv_ms = 5000;
timeout_accept_ms = 5000;
use_hlpscope = true;
use_hlpgetresult = false;

if ischar(func)
    func = str2func(func); end;

% transfer the current scope (optional)
if use_hlpscope
    scope = hlp_resolveall;
    varargin = [{scope,func} varargin];
    func = @hlp_scope;
end

% construct message
fprintf('serializing task data...');
messagebody = hlp_serialize([{1},{func},varargin]);
disp('done.');

% --- set up TCP server to receive the reply ---

fprintf('opening TCP server to receive replies...');
while true
    try
        serv = ServerSocket(2000+round(rand*60000),receive_backlog);
        break;
    catch,end
end
serv.setReceiveBufferSize(receive_buffer_size);
serv.setSoTimeout(timeout_accept_ms);
close_server = onCleanup(@()serv.close());

% determine server's address for replies
return_port = num2str(serv.getLocalPort());
return_host = hlp_hostip;
return_address = [return_host ':' num2str(return_port)];
disp(['done; will receive on ' return_address]);

% --- send task ---

% determine target address
idx = find(address==':',1);
host = address(1:idx-1);
port = str2num(address(idx+1:end));
destination = InetSocketAddress(host,port);

% connect
fprintf('connecting to %s:%i...',host,port);
outconn = Socket();
outconn.connect(destination,timeout_connect_ms); % 3-second timeout
outconn.setTcpNoDelay(1);
outconn.setSoTimeout(timeout_send_ms);
outconn.setSoLinger(true,3);
close_outconn = onCleanup(@()outconn.close());
disp('done.');

% send message
fprintf('now sending (%.0fkb)...',length(messagebody)/1024);
out = DataOutputStream(outconn.getOutputStream());
out.writeInt(taskid);
out.writeInt(length(return_address));
out.writeBytes(return_address);
out.writeInt(length(messagebody));
out.write(messagebody);
out.flush();

% wait for handshake
fprintf('done; now waiting for checkum...');
in = DataInputStream(outconn.getInputStream());
checksum = in.readInt();
fprintf('done; checking...');
if checksum == taskid+length(return_address)+length(messagebody)
    fprintf('all clear');
else
    fprintf('ERROR: unexpected checksum');
end
% close socket
fprintf('; now closing...');
outconn.close();
disp('done.');

% --- receive reply ---

% wait for and get an incoming connection
fprintf('now waiting for reply...');
while true
    try
        conn = serv.accept();
    catch
        continue;
    end
    conn.setSoTimeout(timeout_recv_ms);
    conn.setTcpNoDelay(1);

    % receive reply
    fprintf('done; receiving data.\n');
    in = DataInputStream(conn.getInputStream());
    cr = ChunkReader(in);
    received_taskid = in.readInt();
    if received_taskid ~= taskid
        fprintf('got unexpected task id (%i instead of %i); back to waiting...\n',received_taskid,taskid);
        conn.close();
    else
        disp(['got correct task id: ' num2str(received_taskid)]);
    end
    fprintf('receiving reply...');
    replybody = typecast(cr.readFully(in.readInt()),'uint8');
    fprintf('got reply (%.0fkb).\n',length(replybody)/1024);

    fprintf('now closing...');
    conn.close();
    disp('done.');

    % construct result
    fprintf('parsing result...');
    reply = hlp_deserialize(replybody);
    disp('done.');

    % post-process result
    if isfield(reply{2},{'message','identifier','stack'})
        disp('Got an exception; traceback: ');
        rethrow(reply{2});
    else
        result = reply{2};
    end
    break;
end

if use_hlpgetresult
    varargout = result;
else
    varargout = {result};
end

repo: REST


text:
function score = SSEQ( imdist )

    feature=feature_extract(imdist,3);
    score=SSQA_by_f(feature);

end

author: thomas-koehler

repo: SupER


text:
function [b, mpe] = backT_mpe(engine, f, ev2, t)

bnet = bnet_from_engine(engine);
ss = bnet.nnodes_per_slice;

if t==1
  % ev2 is just the evidence on slice 1
  [mpe, b.clpot] = find_max_config(engine.jtree_engine1, f.clpot, f.seppot, ev2);
else
  [mpe, b.clpot] = find_max_config(engine.jtree_engine, f.clpot, f.seppot, ev2);
  mpe = mpe((1:ss)+ss); % extract values for slice 2
end
for c=1:length(b.clpot)
  [b.clpot{c}, ll(c)] = normalize_pot(b.clpot{c});
end
b.t = t;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% This is an example of how to write the config file
% Line 1 - name of the strategy
% Line 2 - number of Parameters
% Line 3 - Name of Parameter 1
% Line 4 - Default Value of Parameter 1
% ...(do not leave an empty line)
Template Algorithm
3
p1
0
p2
1
p3
10

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function F = tempCompCurve(x,Tdata)
%% Calculate Temperature Curve given Resistor and Thermistor Values
% Copyright (c) 2012, MathWorks, Inc.
%% Input voltage
Vin = 1.1;

%% Thermistor Calculations
% Values in x: R1 R2 R3 R4 RTH1(T_25degc) Beta1 RTH2(T_25degc) Beta2
% Thermistors are represented by:
%   Room temperature is 25degc: T_25
%   Standard value is at 25degc: RTHx_25
%   RTHx is the thermistor resistance at various temperature
% RTH(T) = RTH(T_25degc) / exp (Beta * (T-T_25)/(T*T_25)))
T_25 = 298.15;
T_off = 273.15;
Beta1 = x(6);
Beta2 = x(8);
RTH1 = x(5) ./ exp(Beta1 * ((Tdata+T_off)-T_25)./((Tdata+T_off)*T_25));
RTH2 = x(7) ./ exp(Beta2 * ((Tdata+T_off)-T_25)./((Tdata+T_off)*T_25));

%% Define equivalent circuits for parallel R's and RTH's
R1_eq = x(1)*RTH1./(x(1)+RTH1);
R3_eq = x(3)*RTH2./(x(3)+RTH2);

%% Calculate voltages at Point B
F = Vin * (R3_eq + x(4))./(R1_eq + x(2) + R3_eq + x(4));


text:
function ChangeBinaryGain(filename,varargin)

%ChangeBinaryGain - Change gain in a multiplexed binary file.
%
%  USAGE
%
%    ChangeBinaryGain(filename,<options>)
%
%    filename       file to read
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'output'      output file name (default = append '-gains')
%     'nChannels'   number of data channels in the file (default = 1)
%     'gains'       list of gains for each channel, or one value for all
%     'precision'   sample precision (default = 'int16')
%     'skip'        number of bytes to skip after each value is read
%                   (default = 0)
%    =========================================================================

% Copyright (C) 2004-2011 by Michaël Zugaro
%
% 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.

% Default values
nChannels = 1;
precision = 'int16';
warning(['this function assumes int16 precision, if your file is not int16 use load/resample.m'])
skip = 0;
gains = [];
output = [filename '-gains'];

if nargin < 1 | mod(length(varargin),2) ~= 0,
  error('Incorrect number of parameters (type ''help ChangeBinaryGain'' for details).');
end

% Parse options
for i = 1:2:length(varargin),
	if ~ischar(varargin{i}),
		error(['Parameter ' num2str(i+3) ' is not a property (type ''help ChangeBinaryGain'' for details).']);
	end
	switch(lower(varargin{i})),
		case 'output',
			output = varargin{i+1};
			if ~isstring_FMAT(output),
				error('Incorrect value for property ''output'' (type ''help ChangeBinaryGain'' for details).');
			end
		case 'nchannels',
			nChannels = varargin{i+1};
			if ~isiscalar(nChannels,'>0'),
				error('Incorrect value for property ''nChannels'' (type ''help ChangeBinaryGain'' for details).');
			end
		case 'gains',
			gains = varargin{i+1};
			if ~isdvector(gains,'>=0'),
				error('Incorrect value for property ''gains'' (type ''help ChangeBinaryGain'' for details).');
			end
		case 'precision',
			precision = varargin{i+1};
			if ~isstring_FMAT(precision),
				error('Incorrect value for property ''precision'' (type ''help ChangeBinaryGain'' for details).');
			end
		case 'skip',
			skip = varargin{i+1};
			if ~isiscalar(skip,'>=0'),
				error('Incorrect value for property ''skip'' (type ''help ChangeBinaryGain'' for details).');
			end
		otherwise,
			error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help ChangeBinaryGain'' for details).']);
	end
end

% Check output file
f = fopen(output,'r');
if f ~= -1,
	error(['Output file ''' output ''' exists.']);
end

% Check gains
gains = gains(:)';
nGains = length(gains);
if nGains == 1,
	gains = repmat(gains,1,nChannels);
	nGains = nChannels;
elseif nGains ~= nChannels,
	error('Incorrect number of gains (type ''help ChangeBinaryGain'' for details).');
end

% Read in chunks, apply gains and save
offset = 0;
nSamples = floor(100000/nChannels);
while true,
	data = bz_LoadBinary(filename,'nChannels',nChannels,'offset',offset,'samples',nSamples,'precision',precision,'skip',skip);
	data = data .* repmat(gains,size(data,1),1);
	SaveBinary(output,data,'precision',precision,'mode','append');
	offset = offset + nSamples;
	if size(data,1) < nSamples, break; end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function y = blkdiag( varargin )

%Disciplined convex/geometric programming information for BLKDIAG:
%   BLKDIAG imposes no convexity restrictions on its arguments.

nv = 0;
nz = 0;
sz = [ 0, 0 ];
for k = 1 : nargin,
    x  = cvx( varargin{k} );
    sx = x.size_;
    if length( sx ) > 2,
        error( 'N-D matrices not supported.' );
    end
    b  = x.basis_;
    sz = sz + sx;
    nv = max( nv, size( b, 1 ) );
    nz = nz + nnz( b );
    varargin{k} = x;
end
bz = sparse( [], [], [], prod( sz ), nz, nv );
roff = 0;
coff = 0;
for k = 1 : nargin,
    x  = varargin{k};
    b  = x.basis_;
    sx = x.size_;
    ndxr = ( roff : roff + sx( 1 ) - 1 )';
    ndxr = ndxr( :, ones( 1, sx( 2 ) ) );
    ndxc = ( coff : coff + sx( 2 ) - 1 );
    ndxc = ndxc( ones( 1, sx( 1 ) ), : );
    bz( 1 : size( b, 1 ), ndxc( : ) * sz( 1 ) + ndxr( : ) + 1 ) = b; %#ok
    roff = roff + sx( 1 );
    coff = coff + sx( 2 );
end
y = cvx( sz, bz );

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function disp(F)
%DISP   Display a SPHEREFUNV.
%
% See also SPHEREFUNV/DISPLAY.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

loose = strcmp(get(0,'FormatSpacing'),'loose');

% Compact version:
if ( isempty(F) )
    fprintf('empty spherefunv\n\n')
    return
end

if ( F.isTransposed )
    tString = 'Row vector';
else
    tString = 'Column vector';
end

disp(['   spherefunv object ' '(' tString ') containing' ])
if ( loose )
    fprintf('\n');
end

% Display its two SPHERFUN halves.
for j = 1:3
    disp( F.components{j} );
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function struct = newCERRStructure(scanSet, planC, colorNum)
%"newCERRStructure"
%   Returns an empty, but valid CERR structure. This includes a null string
%   for the structure name and an appropriately formatted contour field.
%
%   scanSet is the index of the scan to register the new structure to.
%   This cannot be changed later.
%
%   Since structure contours require an item for each slice, a planC with 1+ valid
%   scanSets must be in memory.
%   in memory.
%
%JRA 6/4/04
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.
%
%Usage:
%   function struct = newCERRStructure(scanSet, planC)

global stateS

if ~exist('planC', 'var') %wy
    global planC
end

indexS = planC{end};

%Get the structure template from initializeCERR.
struct = initializeCERR('structures');

%Create and populate empty contour field.
try
    nSlices = size(getScanArray(planC{indexS.scan}(scanSet)), 3);
catch
    error('Cannot create new structure if scan is not valid.');
    return
end
segments.points = [];
[contour(1:nSlices).segments] = deal(segments);
[contour(1:nSlices).referencedSopInstanceUID] = deal({planC{indexS.scan}(1).scanInfo(:).sopInstanceUID});
[contour(1:nSlices).referencedSopClassUID] = deal({planC{indexS.scan}(1).scanInfo(:).sopClassUID});

%Create empty name field.
structureName = '';

%Visibility defaults on.
visible = 1;

%Assign color for this structure
assocScanV = getStructureAssociatedScan(1:length(planC{indexS.structures}),planC);
scanIndV = find(assocScanV==scanSet);
if ~exist('colorNum','var')
    colorNum = length(scanIndV) + 1;
end
if ~isfield(stateS,'optS')
%     color = stateS.optS.colorOrder( mod(colorNum-1, size(stateS.optS.colorOrder,1))+1,:);
% else
    stateS.optS = opts4Exe([getCERRPath,'CERROptions.json']);
%     color = stateS.optS.colorOrder( mod(colorNum-1, size(stateS.optS.colorOrder,1))+1,:);
end
colorArr = stateS.optS.colorOrder;
color = setStructureColor(planC,colorArr);
%Assign these values to structure.
struct(1).contour        = contour;
struct(1).structureName  = structureName;
struct(1).visible        = visible;
struct(1).associatedScan = scanSet;
struct(1).assocScanUID   = planC{indexS.scan}(scanSet).scanUID;
struct(1).strUID         = createUID('structure');
struct(1).structureColor = color;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function res = difonzo(Y,x,z,ta,sc,type,f)
% PURPOSE: Multivariate temporal disaggregation with transversal constraint
% ----------------------------------------------------------------------------------
% SYNTAX: res = difonzo(Y,x,z,ta,sc,type,f);
% ----------------------------------------------------------------------------------
% OUTPUT: res: a structure
%         res.meth  = 'Multivariate Di Fonzo';
%         res.N     = Number of low frequency data
%         res.n     = Number of high frequency data
%         res.pred  = Number of extrapolations
%         res.ta    = Type of disaggregation
%         res.sc     = Frequency conversion
%         res.type  = Model for high frequency innovations
%         res.beta  = Model parameters
%         res.y     = High frequency estimate
%         res.d_y   = High frequency estimate: std. deviation
%         res.z     = High frequency constraint
%         res.et    = Elapsed time
% ----------------------------------------------------------------------------------
% INPUT: Y: NxM  ---> M series of low frequency data with N observations
%        x: nxm  ---> m series of high frequency data with n observations, m>=M see (*)
%        z: nzx1 ---> high frequency transversal constraint with nz obs.
%        ta: type of disaggregation
%            ta=1 ---> sum (flow)
%            ta=2 ---> average (index)
%            ta=3 ---> last element (stock) ---> interpolation
%            ta=4 ---> first element (stock) ---> interpolation
%        sc: number of high frequency data points for each low frequency data points
%            sc= 4 ---> annual to quarterly
%            sc=12 ---> annual to monthly
%            sc= 3 ---> quarterly to monthly
%        type: model for the high frequency innvations
%            type=0 ---> multivariate white noise
%            type=1 ---> multivariate random walk
% (*) Optional:
%        f: 1xM ---> Set the number of high frequency indicators linked to
%                    each low frequency variable. If f is explicitly included,
%                    the high frequency indicators should be placed in
%                    consecutive columns
% ----------------------------------------------------------------------------------
% NOTE: Extrapolation is automatically performed when n>sN.
%       If n=nz>sN restricted extrapolation is applied.
%       Finally, if n>nz>sN extrapolation is perfomed in constrained
%       form in the first nz-sN observatons and in free form in
%       the last n-nz observations.
% ----------------------------------------------------------------------------------
% LIBRARY: aggreg, dif, vec, desvec
% ----------------------------------------------------------------------------------
% SEE ALSO: denton, rossi, mtd_print, mtd_plot
% ----------------------------------------------------------------------------------
% REFERENCE: Di Fonzo, T.(1990)"The estimation of M disaggregate time
% series when contemporaneous and temporal aggregates are known", Review
% of Economics and Statistics, vol. 72, n. 1, p. 178-182.

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

% Version 3.0 [August 2006]

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

t0 = clock;

%--------------------------------------------------------
%       Preliminary checking

[N,M] = size(Y);
[n,m] = size(x);
[nz,mz] = size(z);

if ((M > m) | (n < sc*N) | (mz ~= 1) | (nz > n) | (nz < sc*N))
   error (' *** INCORRECT DIMENSIONS *** ');
else
   % Number of extrapolations
   h1 = n - nz;
   h2 = n - sc*N;
   clear nz mz;
end

%--------------------------------------------------------
%       Checking of "ta"

if (ta < 1) | (ta > 4)
    error (' *** INCORRECT TA OPTION *** ');
end

%--------------------------------------------------------
%       Checking of "sc"

if (sc ~= 3) & (sc ~= 4) & (sc ~= 12)
    error (' *** INCORRECT FREQUENCY CONVERSION (sc) *** ');
end

%--------------------------------------------------------
%       Checking of "type"

if (type < 0) | (type > 1)
    error (' *** INCORRECT TYPE OPTION *** ');
end

%--------------------------------------------------------
%       Checking (and definition) of vector f

if (nargin == 6)
    f = ones(1,M);
end

if ( (f < 1) | (sum(f) ~= m) | (length(f) ~= M))
   error (' *** IMPROPER ASSIGNMENT OF INDICATORS IN f VECTOR *** ');
end

%--------------------------------------------------------
%  **** CONSTRAINT MATRICES ***
%--------------------------------------------------------
% Required:
%              H1 ---> transversal
%              H2 ---> longitudinal
%
%---------------------------------------------------------------
%       Generate H1: (n-h1) x nM

H1 = kron(ones(1,M),[eye(n-h1) zeros(n-h1,h1)]);

%---------------------------------------------------------------
%       Generate H2: NM x nM.
%
% Generation of aggregation matrix C

C = aggreg(ta,N,sc);
C = [C zeros(N,h2)];

H2 = kron(eye(M),C);

%---------------------------------------------------------------
%       Generate H: (n-h1+NM) x nM.
%
%       H = [ H1
%             H2 ]

H = [H1
   H2];

%--------------------------------------------------------
%  **** PREPARING DATA MATRICES ***
%--------------------------------------------------------
% Required:
%               x_diag
%               Y_big,  Y_e
%               X_diag, X_e

%--------------------------------------------------------
%       Generate x_diag: nM x M+m
%
% It is a diagonal matrix formed by the high frequency
% indicators, including a vector of ones for the intercept
%
%       x_diag = [ x1 0  0  ... 0
%                  0  x2 0  ... 0
%                  0  0  x3 ... 0
%                  ..............
%                  0  0  0  ... xM ]
%
% It is made by means of a recursion.

ac(1)=f(1);                             % Initialization of the recursion
x_diag = [ones(n,1) x(:,1:ac(1))];

j=2;
while (j <= M)
   xaux = [ones(n,1) x(:,ac(j-1)+1:ac(j-1)+f(j))];
   [a2,b2] = size(xaux);
   [a1,b1] = size(x_diag);
   x_diag = [ x_diag         zeros(a1,b2)
              zeros(a2,b1)   xaux ];
   ac(j) = ac(j-1) + f(j);
   j = j + 1;
end
clear xaux;

%--------------------------------------------------------
%       Generate X_diag: NM x M+m
%
% Low frequency analog of x_diag. It is the result of
% applying the temporal aggregator H2 to x_diag.

X_diag = H2 * x_diag;

%--------------------------------------------------------
%       Generate X_e: (n-h1+NM) x M+m
%
%
% It is the result of applying the complete aggregator H
% (temporal as well as transversal).
% Lower part of X_e is X_diag.

X_e = H * x_diag;

%--------------------------------------------------------
%       Generate Y_big: NM x 1
%
% It is column vector containing all the observations on the
% low frequency series according to: Y_big = [Y1 Y2 ... YM]'
% Formally: Y_big = vec(Y)

Y_big = vec(Y);

%--------------------------------------------------------
%       Generate Y_e: (n-h1+NM) x 1
%
% It is column vector containing the transversal constraint
% and all the observations on the low frequency series
% according to: Y_e = [ z Y1 Y2 ... YM]' = [z Y_big]'

Y_e = [ z
      Y_big];

%--------------------------------------------------------
%  **** PRELIMINARY ESTIMATION OF SIGMA ***
%--------------------------------------------------------
% The method of di Fonzo requires the previous estimation of VCV
% matrix SIGMA for the (implied) low frequency model. This
% preliminary estimation is performed by means of estimating, equation
% by equation, the model. Formally, this is equivalent to estimate an
% unrelated SURE model. Computationally, this is also the applied procedure.

BETA = (X_diag' * X_diag) \ (X_diag' * Y_big); % OLS estimator
U_big = Y_big - X_diag * BETA;                 % Residuals in vec format

% Residuals (columnwise) U: NxM

U = desvec(U_big,M);

% Preliminary estimation of SIGMA

SIGMA = cov(U,1);

%--------------------------------------------------------
%  **** APPLYING DI FONZO PROCEDURE ***
%--------------------------------------------------------

%--------------------------------------------------------
%       High frequency VCV matrix v: nM x nM

switch type
case 0 % White noise
   v = kron(SIGMA,eye(n));
case 1 % Random walk, with U(0)=0
   D = dif(1,n);
   DDi = inv(D'*D);
   v = kron(SIGMA,DDi);
end;

%--------------------------------------------------------
%       Low frequency VCV matrix V: (n-h1+NM) x (n-h1+NM)
%       and its generalized inverse

V = H * v * H';
Vi = pinv(V);      % Moore-Penrose generalized inverse

%--------------------------------------------------------
%       Generation of distribution filter L: nM x (n-h1+NM)

L = v * H' * Vi;

%--------------------------------------------------------
%	     GLS estimation of beta in a SURE context

beta = (X_e' * Vi * X_e) \ (X_e' * Vi * Y_e);

U_e = Y_e - X_e * beta;

%--------------------------------------------------------
%       Estimation of high frequency series

y_big = x_diag * beta + L * U_e;

% Series y columnwise y: nxM

y = desvec(y_big,M);

%--------------------------------------------------------
%       VCV matrix of estimations y: nM x nM

sigma_y = (eye(n*M) - L*H)*v + ...
   (x_diag - L*X_e)*inv(X_e'*Vi*X_e)*(x_diag - L*X_e)';

% Vector format of std. dev.

d_y_big = sqrt(diag(sigma_y));

% Std. dev. series in column format dt_y: n x M

d_y = desvec(d_y_big,M);

% -----------------------------------------------------------------------
% Loading the structure
% -----------------------------------------------------------------------
% Basic parameters

res.meth = 'Multivariate Di Fonzo';
res.N = N;
res.n = n;
res.pred = h2;
res.ta= ta;
res.sc = sc;
res.type = type;

% -----------------------------------------------------------------------
% Parameters

res.beta=beta;

% -----------------------------------------------------------------------
% Series

res.y   = y;
res.d_y = d_y;
res.z   = z;

% -----------------------------------------------------------------------
% Elapsed time

res.et        = etime(clock,t0);


text:
function article=indefiniteArticle(phrase)
% A simple function that returns the indefinite articles "a" or "an" based on a given word or phrase.
% based on the javascript implementation: https://github.com/rigoneri/indefinite-article.js/blob/master/indefinite-article.js
% indefiniteArticle('')
% indefiniteArticle('apple pie')
% indefiniteArticle('vision group')
% indefiniteArticle('Priceton local restaurant')

   % Getting the first word
    [~,match] = regexp(phrase,'\w+');

    if ~isempty(match)
        word = phrase(1:match(1));
    else
        article=''; return;
    end

    l_word = lower(word);
    % Specific start of words that should be preceeded by 'an'
    alt_cases = {'honest', 'hour', 'hono'};
    for i=1:length(alt_cases)
        t = strfind(l_word, alt_cases{i});
        if (~isempty(t) && t(1) == 1)
            article='an'; return;
        end
    end

    % Single letter word which should be preceeded by 'an'
    if length(l_word) == 1
        if ~isempty(strfind('aedhilmnorsx',l_word))
            article='an'; return;
        else
            article='a'; return;
        end
    end

    % Capital words which should likely be preceeded by 'an'

    match = regexp(word,'(?!FJO|[HLMNS]Y.|RY[EO]|SQU|(F[LR]?|[HL]|MN?|N|RH?|S[CHKLMNPTVW]?|X(YL)?)[AEIOU])[FHLMNRSX][A-Z]', 'once');
    if ~isempty(match)
        article='an'; return;
    end

    % Special cases where a word that begins with a vowel should be preceeded by 'a'
    regexes = {'^e[uw]', '^onc?e\b', '^uni([^nmd]|mo)', '^u[bcfhjkqrst][aeiou]'};
    for i=1:length(regexes)
        match = regexp(l_word,regexes{i}, 'once');
        if ~isempty(match)
            article='a'; return;
        end
    end

    % Special capital words (UK, UN)
    match = regexp(word,'^U[NK][AIEO]', 'once');
    if ~isempty(match)
        article='a'; return;
    else
        if strcmp(word, upper(word))
            if ~isempty(strfind('aedhilmnorsx',l_word(1)))
                article='an'; return;
            else
                article='a'; return;
            end
        end
    end

    % Basic method of words that begin with a vowel being preceeded by 'an'
    if ~isempty(strfind('aeiou',l_word(1)))
        article='an'; return;
    end

    % Instances where y follwed by specific letters is preceeded by 'an'
    match = regexp(l_word,'^y(b[lor]|cl[ea]|fere|gg|p[ios]|rou|tt)', 'once');
    if ~isempty(match)
        article='an'; return;
    end

    article='a'; return;

author: jianxiongxiao

repo: ProfXkit


text:
% hungarianLSAPE.m Help file for hungarianLSAP MEX-file.
%  hungarianLSAPE.cpp - Compute a solution to (symmetric or asymmetric) LSAP with Hungarian algorithm
%
%      [rho,varrho] = hungarianLSAPE(C,init_type,forb)
%         [rho,u,v] = hungarianLSAPE(C,init_type,forb)
%  [rho,varrho,u,v] = hungarianLSAPE(C,init_type,forb)
%
%  Given a nxm cost matrix C (integer ou floatting values)
%  with last column encoding removal costs and last row insertion costs, it computes:
%
%  - a solution rho to the LSAPE, i.e. a mapping from the rows of C to its columns,
%    and the mapping varrho from the columns to the rows
%
%    rho is represented by a (n-1)x1 matrix so that rho(i)=j means that i is assigned to j
%    varrho is a 1x(m-1) matrix so that varrho(j)=i means that j is assigned to i
%    rho(i)=m or varrho(j)=n encode assignments to the null element (removal or insertion)
%
%  - a solution (u,v) to its dual problem (labeling problem)
%
%  optional init_type:
%    0: no initialization (u=v=0)
%    1 (default): classical initialization (u=min(C) on rows, v=min(C-u) on columns)
%
%  optional boolean parameter forb:
%    true  -> forbidden assignments are represented by negative cost values
%    false -> no forbidden assignments (by default)
%
%  This is a MEX-file for MATLAB.
%  This file is part of LSAPE.
%  LSAPE is free software: you can redistribute it and/or modify it
%  under the terms of the CeCILL-C License. See README for more details.
%
%     Copyright 2015-2017
%      authors: Sebastien Bougleux
%  institution: Normandie Univ, CNRS - ENSICAEN - UNICAEN, GREYC, France
%   last modif: July 5 2017
%

author: optimaltransport


text:
 function y = fatrix2_def_forw(M, x)
%function y = fatrix2_def_forw(M, x)
%| basic default 'forw' function: M * x
y = M * x;

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function varargout = grMergeMultipleEdges(nodes, edges)
%GRMERGEMULTIPLEEDGES Remove all edges sharing the same extremities
%
%   [NODES2, EDGES2] = grMergeMultipleEdges(NODES, EDGES)
%   Remove configuration with two edges sharing the same 2 nodes.
%

%   -----
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 13/08/2003.
%

%   HISTORY
%   10/02/2004 doc
%   2011-05-18 rename to grMergeMultipleEdges

rmedge = [];
for e = 1:length(edges)
    edge = edges(e, :);
    for e2 = e+1:length(edges)
        if (edge(1) == edges(e2, 1) && edge(2) == edges(e2, 2)) || ...
           (edge(1) == edges(e2, 2) && edge(2) == edges(e2, 1))
                rmedge(length(rmedge)+1) = e2; %#ok<AGROW>
        end
    end
end

[nodes, edges] = grRemoveEdges(nodes, edges, rmedge);

% process output depending on how many arguments are needed
if nargout == 1
    out{1} = nodes;
    out{2} = edges;
    varargout{1} = out;
end

if nargout == 2
    varargout = {nodes, edges};
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function make_fig2(A, y, cuckoo_idx, splunk_idx, t_created, virus, names)

fold = 2;
cutoff = .3;

I_1 = cv_part(y, fold, cutoff, 'standard', t_created, cuckoo_idx, splunk_idx);
I_2 = cv_part(y, fold, cutoff, 'creation', t_created, cuckoo_idx, splunk_idx);
I_3 = cv_part(y, fold, cutoff, 'type', t_created, cuckoo_idx, splunk_idx, virus);

h = figure;
set(h,'Renderer','painters','Position',[100 900 1000 400]);

subaxis(2,4,1,'Margin',0.00,'Padding',-0.00,'PaddingBottom',0.0,'PaddingTop',0.02,'PaddingLeft',0.02,'PaddingRight',-0.01,'MarginBottom',0.08,'MarginLeft',0.03,'MarginRight',0.10,'MarginTop',0.03);

ha = subaxis(2,4,1,1);
[pred02, pred05, logit_model1] = run(A, y, I_1, false, false, t_created, cuckoo_idx, splunk_idx);
pos1 = get(gca,'Position');
annotation(h,'textbox',[pos1(1)-0.05,pos1(2)+pos1(4)-0.01,.1,.05],'String','(A)','fontsize',14,'EdgeColor','none');
ylabel('TPR');

ha = subaxis(2,4,2,1);
[~, ~, logit_model2] = run(A, y, I_2, false, true, t_created, cuckoo_idx, splunk_idx);
pos1 = get(gca,'Position');
annotation(h,'textbox',[pos1(1)-0.05,pos1(2)+pos1(4)-0.01,.1,.05],'String','(B)','fontsize',14,'EdgeColor','none');

ha = subaxis(2,4,3,1);
[~, ~, logit_model3] = run(A, y, I_3, false, false, t_created, cuckoo_idx, splunk_idx);
pos1 = get(gca,'Position');
annotation(h,'textbox',[pos1(1)-0.05,pos1(2)+pos1(4)-0.01,.1,.05],'String','(C)','fontsize',14,'EdgeColor','none');

ha = subaxis(2,4,4,1);
[h1, h2] = virus_names(y, pred02, pred05, virus);
pos1 = get(gca,'Position');
annotation(h,'textbox',[pos1(1)-0.05,pos1(2)+pos1(4)-0.01,.1,.05],'String','(D)','fontsize',14,'EdgeColor','none');
pos1 = get(h1,'Position');
set(h1,'Position', [pos1(1)+0.05, pos1(2), pos1(3), pos1(4)]);
if (h1~=h2)
    pos1 = get(h2,'Position');
    set(h2,'Position', [pos1(1)+0.05, pos1(2), pos1(3), pos1(4)]);

end

ha = subaxis(2,4,1,2);
[pred02, pred05] = run(A, y, I_1, true, false, t_created, cuckoo_idx, splunk_idx, logit_model1);
pos1 = get(gca,'Position');
annotation(h,'textbox',[pos1(1)-0.05,pos1(2)+pos1(4)-0.01,.1,.05],'String','(E)','fontsize',14,'EdgeColor','none');
ylabel('TPR');
xlabel('FPR');

ha = subaxis(2,4,2,2);
run(A, y, I_2, true, true, t_created, cuckoo_idx, splunk_idx, logit_model2);
pos1 = get(gca,'Position');
annotation(h,'textbox',[pos1(1)-0.05,pos1(2)+pos1(4)-0.01,.1,.05],'String','(F)','fontsize',14,'EdgeColor','none');
xlabel('FPR');

ha = subaxis(2,4,3,2);
run(A, y, I_3, true, false, t_created, cuckoo_idx, splunk_idx, logit_model3);
pos1 = get(gca,'Position');
annotation(h,'textbox',[pos1(1)-0.05,pos1(2)+pos1(4)-0.01,.1,.05],'String','(G)','fontsize',14,'EdgeColor','none');
xlabel('FPR');

ha = subaxis(2,4,4,2);
[h1, h2] = virus_names(y, pred02, pred05, virus);
pos1 = get(gca,'Position');
annotation(h,'textbox',[pos1(1)-0.05,pos1(2)+pos1(4)-0.01,.1,.05],'String','(H)','fontsize',14,'EdgeColor','none');
pos1 = get(h1,'Position');
set(h1,'Position', [pos1(1)+0.05, pos1(2), pos1(3), pos1(4)]);
pos1 = get(h2,'Position');
if (h1~=h2)
    set(h2,'Position', [pos1(1)+0.05, pos1(2), pos1(3), pos1(4)]);
end
xlabel('Fraction Detected');

end

function Ic = valid_cuckoo_only(I, cuckoo_idx)

    Ic = {};
    for iter=1:length(I)
      Ic{iter}.training = I{iter}.training | (I{iter}.test & ~cuckoo_idx);
      Ic{iter}.test = I{iter}.test & cuckoo_idx;
   end
end

function [pred02, pred05, logit_model_out] = run(A, y, I, cuckoo_only, time_split, t_created, cuckoo_idx, splunk_idx, logit_model_in)

    A = double(A>0);
    y = double(y>0);

    logit = init();
    bagged = init();
    logit2 = init();
    bagged2 = init();
    logit3 = init();
    bagged3 = init();

    pred02 = NaN*ones(size(y));
    pred05 = NaN*ones(size(y));

    %compute valid
    valid = false(size(y));
    for iter=1:length(I)
        valid = valid | I{iter}.training | I{iter}.test;
    end

    logit_model_out = {};

    for iter=1:length(I)

        I_train = I{iter}.training;
        I_test = I{iter}.test;

        if (time_split)

            t_smallest = min(t_created(I_test & t_created>0 & y>0 & cuckoo_idx));

            I_train = valid & I_train & (y==0 | ~cuckoo_idx | t_created<(t_smallest-86400 * 365 *2));
            I_test = valid & ~I_train & (y==0 | ~cuckoo_idx | t_created<(t_smallest-86400 * 365 *2));
        end

        if (sum(I_train & I_test)>0)
            error('Testing and training intersecting.');
        end

        A_train = A(I_train, :);
        y_train = y(I_train, :);

        if (nargin<9)
            [logit_model, boost_model] = run_pipeline(A_train, y_train);
            logit_model_out{iter} = logit_model;
        else
            logit_model = logit_model_in{iter};
            boost_model = logit_model;
            logit_model_out{iter} = logit_model;
        end

        %graft data
        if (cuckoo_only)
            [A_test, y_test, I_test] = gen_graft_data(A,y,cuckoo_idx,splunk_idx, I_test, logit_model);
        else
            I_test = I_test & ~splunk_idx;
            A_test = A(I_test, :);
            y_test = y(I_test, :);
        end

        [logit, T02, T05] = record_run(A_test, y_test, logit_model, logit);
        [bagged] = record_run(A_test, y_test, boost_model, bagged);

        %record the predictions
        pred02(I_test) = logit_model.predict(A_test)>=T02;
        pred05(I_test) = logit_model.predict(A_test)>=T05;

        if (time_split)

            I_test = valid & ~I_train & (y==0 | ~cuckoo_idx | t_created>=(t_smallest-86400 * 365 * 1));

           if (sum(I_train & I_test)>0)
                error('Testing and training intersecting.');
            end

            %A_train = A(I_train, :);
            %y_train = y(I_train, :);

            %[logit_model2, boost_model2] = run_pipeline(A_train, y_train);
            logit_model2 = logit_model;
            boost_model2 = boost_model;

            %graft data
            if (cuckoo_only)
                [A_test, y_test,I_test] = gen_graft_data(A,y,cuckoo_idx,splunk_idx, I_test, logit_model2);
            else
                I_test = I_test & ~splunk_idx;
                A_test = A(I_test, :);
                y_test = y(I_test, :);
            end

            [logit2] = record_run(A_test, y_test, logit_model2, logit2);
            [bagged2] = record_run(A_test, y_test, boost_model2, bagged2);

            I_test = valid & ~I_train & (y==0 | ~cuckoo_idx | t_created>=(t_smallest-86400 * 365 * 0));

            if (sum(I_train & I_test)>0)
                error('Testing and training intersecting.');
            end

            %A_train = A(I_train, :);
            %y_train = y(I_train, :);

            %[logit_model3, boost_model3] = run_pipeline(A_train, y_train);
            logit_model3 = logit_model;
            boost_model3 = boost_model;

            %graft data
            if (cuckoo_only)
                [A_test, y_test,I_test] = gen_graft_data(A,y,cuckoo_idx,splunk_idx, I_test, logit_model3);
            else
                I_test = I_test & ~splunk_idx;
                A_test = A(I_test, :);
                y_test = y(I_test, :);
            end

            [logit3] = record_run(A_test, y_test, logit_model3, logit3);
            [bagged3] = record_run(A_test, y_test, boost_model3, bagged3);

        end

    end

    if (time_split)
        semilogx(mean(logit.X) , mean(logit.Y),'-k', mean(logit2.X) , mean(logit2.Y),'-b', mean(logit3.X) , mean(logit3.Y),'-r', mean(bagged.X) , mean(bagged.Y), '--k', mean(bagged2.X) , mean(bagged2.Y), '--b', mean(bagged3.X) , mean(bagged3.Y), '--r');
        legend(sprintf('0 year, AUC=%.2f', mean(logit.AUC)), sprintf('1 year, AUC=%.2f', mean(logit2.AUC)), sprintf('2 year, AUC=%.2f', mean(logit3.AUC)), 'Location','SouthEast');
        %plot(mean(logit.X) , mean(logit.Y),'-k', mean(logit2.X) , mean(logit2.Y),'-r');
        %legend(sprintf('LR, AUC=%.2f, 1 year', mean(logit.AUC)), sprintf('LR, AUC=%.2f, 2 years', mean(logit2.AUC)), 'Location','SouthEast');
    else

        semilogx(mean(logit.X) , mean(logit.Y),'-k', mean(bagged.X) , mean(bagged.Y), '--k');
        legend(sprintf('AUC=%.2f', mean(logit.AUC)), 'Location','SouthEast');
        %plot(mean(logit.X) , mean(logit.Y),'-k');
    end

    xlim([2e-4, 0.1]);
    ylim([0, 1.0]);
    set(gca, 'XMinorTick','on','YMinorTick','on','TickLength',[0.025 0.025]);

end

function [base] = init()
    base.X = [];
    base.Y = [];
    base.AUC = [];
    base.AUC_splunk = [];
end

function [base,T02, T05] = record_run(A, y, model, base)

    [X,Y,T,auc] = model.getROCPoints(A, y);

    T02 = T(find(X>=.01,1,'first'));
    T05 = T(find(X>=.001,1,'first'));
    AUC_splunk = 0;

    base.X(end+1,:) = X;
    base.Y(end+1,:) = Y;
    base.AUC(end+1) = auc;
    base.AUC_splunk(end+1) = AUC_splunk;

end

 function [h1, h2] = virus_names(y, pred02, pred05, virus)

    if (sum(~isnan(pred02))==0)
        h1 = gca;
        %h2 = axes();
        h2 = gca;
        return;
    end

  for iter=1:length(virus)
      if (strcmp(virus{iter}, 'Trojan.Win32.Generic'))
          virus{iter} = 'Generic';
      else
        x = strsplit(virus{iter},'.');
        virus{iter} = lower(x{1});
        x = strsplit(virus{iter},'-');
        virus{iter} = lower(x{1});
      end
  end

  [unique_names,~,z] = unique(virus);
  d = hist(max(0, z(y>0 & ~isnan(pred02))),length(unique_names));

  d_after = hist(z(y>0 & pred05>0),length(unique_names));
  d_after2 = hist(z(y>0 & pred02>0),length(unique_names));

  [~,I] = sort(d,'descend');

  I_label = [];
  for iter=1:length(I)
    if (strcmp(unique_names(I(iter)),'<nolabel>')==0 && strcmp(unique_names(I(iter)),'Generic')==0)
      I_label = [I_label I(iter)];
    end

    if (length(I_label)>=10)
      break;
    end
  end

  barh([d_after(I_label)'./d(I_label)', d_after2(I_label)'./d(I_label)']);

  hold on;

  set(gca, 'YDir','reverse');
  set(gca,'YTickLabel',unique_names(I_label), 'XMinorTick','on','TickLength',[0.025 0.025]);
  ylabel('Kaspersky Classification');
  xlim([0, 1.0]);
  h1 = gca;

  %legend('FPR=10^{-2}', 'FPR=10^{-3}','Location','SouthWest' );
  %legend('boxoff');

  h2 = axes('Position',get(gca,'Position'),...
              'XAxisLocation','top',...
              'YAxisLocation','right',...
              'Color','none');

  set(h2, 'Ylim', get(h1,'Ylim'), 'YDir','reverse', 'XTickLabel', [], 'YTickLabel', d(I_label), 'YTick', 1:length(I_label));
  ylabel(h2, '# Observations');
  axes(h1);

 end

author: konstantinberlin

repo: malware-windows-audit-log-detection


text:
classdef FollowWall < simiam.controller.Controller

% Copyright (C) 2013, Georgia Tech Research Corporation
% see the LICENSE file included with this software

    properties

        % memory banks
        E_k
        e_k_1

        % gains
        Kp
        Ki
        Kd

        % plot support
        p

        % sensor geometry
        calibrated
        sensor_placement

    end

    properties (Constant)
        inputs = struct('v', 0, 'direction', 'right');
        outputs = struct('v', 0, 'w', 0)
    end

    methods

        function obj = FollowWall()
            obj = obj@simiam.controller.Controller('follow_wall');
            obj.calibrated = false;

            obj.Kp = 2;
            obj.Ki = 0;
            obj.Kd = 0;

            obj.E_k = 0;
            obj.e_k_1 = 0;

%             obj.p = simiam.util.Plotter();
        end

        function outputs = execute(obj, robot, state_estimate, inputs, dt)

            % Compute the placement of the sensors
            if(~obj.calibrated)
                obj.set_sensor_geometry(robot);
            end

            % Unpack state estimate
            [x, y, theta] = state_estimate.unpack();

            % Poll the current IR sensor values 1-5
            ir_distances = robot.get_ir_distances();

            % Interpret the IR sensor measurements geometrically
            ir_distances_wf = obj.apply_sensor_geometry(ir_distances, state_estimate);

            % Compute the heading vector
            d_fw = inputs.d_fw;

            % 1. Select p_2 and p_1, then compute u_fw_t
            if(strcmp(inputs.direction,'right'))
                % Pick two of the right sensors based on ir_distances
                S = [1:3 ; ir_distances(5:-1:3)'];
                [Y,i] = sort(S(2,:));
                S = S(1,i);

                Sp = 5:-1:3;

                S1 = Sp(S(1));
                S2 = Sp(S(2));

                if(S1 < S2)
                    p_1 = ir_distances_wf(:,S2);
                    p_2 = ir_distances_wf(:,S1);
                else
                    p_1 = ir_distances_wf(:,S1);
                    p_2 = ir_distances_wf(:,S2);
                end

            else
                % Pick two of the left sensors based on ir_distances
                S = [1:3 ; ir_distances(1:3)'];
                [Y,i] = sort(S(2,:));
                S = S(1,i);

                if(S(1) > S(2))
                    p_1 = ir_distances_wf(:,S(2));
                    p_2 = ir_distances_wf(:,S(1));
                else
                    p_1 = ir_distances_wf(:,S(1));
                    p_2 = ir_distances_wf(:,S(2));
                end
            end

            u_fw_t = p_2-p_1;

            % 2. Compute u_a, u_p, and u_fw_tp to compute u_fw_p

            u_fw_tp = u_fw_t/norm(u_fw_t);
            u_a = p_1;
            u_p = [x;y];

            u_fw_p = ((u_a-u_p)-((u_a-u_p)'*u_fw_tp)*u_fw_tp);

            % 3. Combine u_fw_tp and u_fw_pp into u_fw;
            u_fw_pp = u_fw_p/norm(u_fw_p);
            u_fw = d_fw*u_fw_tp+(u_fw_p-d_fw*u_fw_pp);

            % Compute the heading and error for the PID controller
            theta_fw = atan2(u_fw(2),u_fw(1));
            e_k = theta_fw-theta;
            e_k = atan2(sin(e_k),cos(e_k));

            e_P = e_k;
            e_I = obj.E_k + e_k*dt;
            e_D = (e_k-obj.e_k_1)/dt;

            % PID control on w
            v = inputs.v;
            w = obj.Kp*e_P + obj.Ki*e_I + obj.Kd*e_D;

            % Save errors for next time step
            obj.E_k = e_I;
            obj.e_k_1 = e_k;

            % plot
%             obj.p.plot_2d_ref(dt, atan2(sin(theta),cos(theta)), theta_fw, 'c');

%             fprintf('(v,w) = (%0.4g,%0.4g)\n', v,w);

            outputs.v = v;
            outputs.w = w;
        end

        % Helper functions

        function ir_distances_wf = apply_sensor_geometry(obj, ir_distances, state_estimate)

            % Apply the transformation to robot frame.

            ir_distances_rf = zeros(3,5);
            for i=1:5
                x_s = obj.sensor_placement(1,i);
                y_s = obj.sensor_placement(2,i);
                theta_s = obj.sensor_placement(3,i);

                R = obj.get_transformation_matrix(x_s,y_s,theta_s);
                ir_distances_rf(:,i) = R*[ir_distances(i); 0; 1];
            end

            % Apply the transformation to world frame.

            [x,y,theta] = state_estimate.unpack();

            R = obj.get_transformation_matrix(x,y,theta);
            ir_distances_wf = R*ir_distances_rf;

            ir_distances_wf = ir_distances_wf(1:2,:);
        end

        function set_sensor_geometry(obj, robot)
            obj.sensor_placement = zeros(3,5);
            for i=1:5
                [x, y, theta] = robot.ir_array(i).location.unpack();
                obj.sensor_placement(:,i) = [x; y; theta];
            end
            obj.calibrated = true;
        end

        function R = get_transformation_matrix(obj, x, y, theta)
            R = [cos(theta) -sin(theta) x; sin(theta) cos(theta) y; 0 0 1];
        end

        function reset(obj)
            % Reset accumulated and previous error
            obj.E_k = 0;
            obj.e_k_1 = 0;
        end

    end

end

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text:
function[image, oriImSize] = e2s2_prepareImage(net, image, maxImageSize)
% [image, oriImSize] = e2s2_prepareImage(net, image, maxImageSize)
%
% Resize the image and subtract the mean image.
%
% Copyright by Holger Caesar, 2015

% Resize image
oriImSize = size(image);
resizeFactor = maxImageSize / max(oriImSize(1:2));
targetSize = ceil(oriImSize(1:2) .* resizeFactor); % ceil corresponds to Matlab's imresize behavior
image = imresize(image, resizeFactor);
assert(size(image, 1) == targetSize(1) && size(image, 2) == targetSize(2));

if numel(net.meta.normalization.averageImage) == 3,
    % Subtract fixed number from each channel
    image(:, :, 1) = image(:, :, 1) - net.meta.normalization.averageImage(1);
    image(:, :, 2) = image(:, :, 2) - net.meta.normalization.averageImage(2);
    image(:, :, 3) = image(:, :, 3) - net.meta.normalization.averageImage(3);
else
    % Resize averageImage and subtract it from each image
    % Note: This cannot be done on the gpu as Matlab's gpu-compatible
    % imresize function can only resize by a constant factor and the image
    % might not be square.
    averageImage = net.meta.normalization.averageImage ./ 255;
    averageImage = imresize(averageImage, targetSize);
    image = image - averageImage;
end;

author: nightrome

repo: matconvnet-calvin


text:
function tran06_values_test ( )

%*****************************************************************************80
%
%% TRAN06_VALUES_TEST demonstrates the use of TRAN06_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRAN06_VALUES_TEST:\n' );
  fprintf ( 1, '  TRAN06_VALUES stores values of\n' );
  fprintf ( 1, '  the transportation function of order 6.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X               FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

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

    if ( n_data == 0 )
      break
    end

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

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function p = predict(Theta1, Theta2, X)
%PREDICT Predict the label of an input given a trained neural network
%   p = PREDICT(Theta1, Theta2, X) outputs the predicted label of X given the
%   trained weights of a neural network (Theta1, Theta2)

% Useful values
m = size(X, 1);
num_labels = size(Theta2, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned neural network. You should set p to a
%               vector containing labels between 1 to num_labels.
%
% Hint: The max function might come in useful. In particular, the max
%       function can also return the index of the max element, for more
%       information see 'help max'. If your examples are in rows, then, you
%       can use max(A, [], 2) to obtain the max for each row.
%

temp1 = [ones(m, 1) X];
temp2 = [ones(m, 1) sigmoid(temp1 * Theta1')];
temp3 = sigmoid(temp2 * Theta2');
[maxTemp3, maxTemp3_2] = max(temp3');
p = maxTemp3_2';

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

end

repo: MachineLearning


text:
function pred = ml_predictqda(trials, model)
% Prediction function for Quadratic Discriminant Analysis.
% Prediction = ml_predictqda(Trials, Model)
%
% In:
%   Trials  : the data a matrix, as in ml_predict
%
%   Model   : predictive model as produced by ml_trainqda
%
% Out:
%   Prediction  : discrete probability distribution, formatted as
%                 {'disc' [NxC] [Cx1]}, with element #2 being the per-class probability and
%                 element #3 the original target values per class
%                 thus, the expected target values are Prediction{2}*Prediction{3}
%
% Examples:
%   targets might look like this: [-1 -1 1 -1 1 -1 -1 1 -1 -1 1 -1 -1 1 ...]'
%
%   model = ml_trainqda(data,targets)
%   p = ml_predictqda(data, model); expectation = p{2}*p{3};
%   now expectation might look like this: [-0.6 -0.9 +0.4 -0.7 +0.8 -0.1 +0.5 +1.0 -0.9 +1.0 -1.0 -1.0 +1.0 ...]'
%
% See also:
%   ml_trainqda
%
%                           Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                           2010-04-03

if isfield(model,'voted')
    % dispatch to the voter
    pred = ml_predictvote(trials,model);
else
    % pre-prune features
    trials = trials(:,model.featuremask);
    % pre-scale the data
    trials = hlp_applyscaling(trials,model.sc_info);
    % calculate the labels
    raw_labels = min(+1,max(-1,-(sum(((trials*model.q).*trials)') + model.l*trials' - model.c)))'; %#ok<UDIM>
    pred = {'disc', [(1-raw_labels)/2 1-(1-raw_labels)/2], model.classes};
end


text:
% MPC parameters
N  = 13;                        % Prediction horizon (number of iterations)
Nu  = N;                        % Control horizon (number of iterations)

Q = [25 0 0];                   % State weights
R = 0.05;                       % du weights, 0
Ru = 0.05;                      % u weights, 0.1
LB = -0.3*ones(Nu,1);           % Lower bound of control input, -0.3
UB = 0.5*ones(Nu,1);            % Upper bound of control input, 0.5
LBdu = nan;                     % Lower bound of control input rate, -0.1
UBdu = nan;                     % Upper bound of control input rate, 0.1
LBo = [-0.2,-1,-1];             % Lower bound of output
UBo = [0.4,1,1];                % Upper bound of output

% Reference trajectory
time_control = 0:dt:Duration;
r = 0.4*(-0.5./(1+exp(time_control./0.1-8)) + 1./(1+exp(time_control./0.1-30)) - 0.4);
figure,plot(time_control,r,'-k'), hold on

xrefFUN = @(t) 0.4*(-0.5./(1+exp(t/0.1-8)) + 1./(1+exp(t/0.1-30)) - 0.4);
plot(time_control,xrefFUN(time_control),'--r')

author: eurika-kaiser

repo: SINDY-MPC


text:
function str=onIff(test)

if test,
  str='on';
else
  str='off';
end

return

author: kristinbranson

repo: JAABA


text:
function [M_Gamma2] = gamma_square (Gamma)

M_Gamma2 = zeros(size(Gamma));
for k=1:size(Gamma,2)
    M_Gamma2(:,k) = Gamma(:,k).^2 + Gamma(:,k).*(1-Gamma(:,k));
end;

author: OHBA-analysis

repo: HMM-MAR


text:
%% For CUDA functionality before running matlab, run
% export LD_LIBRARY_PATH=/pkgs_local/cuda-4.0/lib64/
% export PATH=$PATH:/pkgs_local/cuda-4.0/bin/
% export CUDA_PATH='/pkgs_local/cuda-4.0';

%% If running multiple instances of code on the same machine the number of threads for openmp should be restricted by
% export OMP_NUM_THREADS=4;

%% For Jacket GPU functionality
% addpath ~norouzi/research/jacket/engine;
% ginfo      % shows the gpu information
% gselect(1) % selects the first gpu

% Compilation of the mex files
mex loss_aug_inf_triplet_mex0.cpp CXXFLAGS="\$CXXFLAGS -Wall -fopenmp" LDFLAGS="\$LDFLAGS -fopenmp"
mex loss_aug_inf_mex0.cpp CXXFLAGS="\$CXXFLAGS -Wall -fopenmp" LDFLAGS="\$LDFLAGS -fopenmp"
mex utils/hammknn_mex.cpp utils/linscan.cpp CXXFLAGS="\$CXXFLAGS -fopenmp" LDFLAGS="\$LDFLAGS -fopenmp";

% adding some stuff to the path.
addpath utils; % You should be at thre root of the project to run this line.
addpath /YAEL_HOME/matlab/ % download and install yael library and change this line

% parameters & RUN2
DATA = 'mnist';  % either 'mnist' or 'cifar10'
REG = 1;         % always 1 -- the term that enforces bits to have a mean of zero
NB = 32;         % number of bits
LOSS = 1;        % loss type: 0 -> pairwise hinge / 1 -> triplet ranking
LSCALE = 1;      % An scaling parameter multiplied by the loss
NN = 1000;       % Number of same-class items to be used in the positive set for training.
                 % NN=Inf -> all of the items from the same class
                 % NN=1000 -> find same-class 1000-NN in euclidean distance and set them as target positive examples for each item
GPU = 0;         % GPU=0 -> No GPU / GPU=1 -> GPU (change this to 1 if you have jacket configured)
SHRINK_W = [3e-2 3e-3 3e-4];
                 % weight decay parameter -- cross validating on this parameter on 3 candidates
RUN2             % given the above parameters runs the code

% Linear hash functions on CIFAR-10
DATA = 'cifar10'; REG = 1; NB = 128; LOSS = 1; LSCALE = 1; NN = Inf; GPU = 0; SHRINK_W = 1e-3;
RUN2

% Neural Net because NB is two dimensional. Calles compute_NN_output and compute_NN_grad
DATA = 'mnist'; REG = 1; NB = [512 32]; LOSS = 1; LSCALE = 1; NN = Inf; GPU = 0; SHRINK_W = 1e-4*[1; 1];
RUN2

author: norouzi

repo: hdml

save_path: github-repos/MATLAB/norouzi-hdml


text:
% VL_QUICKSHIFT Quick shift image segmentation
%   Quick shift is a mode seeking algorithm which links each pixel to
%   its nearest neighbor which has an increase in the estimate of the
%   density. These links form a tree, where the root of the tree is
%   the pixel which correspond to the highest mode in the image.
%
%   [MAP,GAPS] = VL_QUICKSHIFT(I, KERNELSIZE, MAXDIST) computes quick shift on the
%   image I. KERNELSIZE is the bandwidth of the Parzen window estimator of
%   the density. Since searching over all pixels for the nearest
%   neighbor which increases the density would be prohibitively
%   expensive, MAXDIST controls the maximum L2 distance between neighbors
%   that should be linked. MAP and GAP represent the resulting forest
%   of trees. They are array of the same size of I.  Each element
%   (pixel) of MAP is and index to the parent elemen in the forest and
%   GAP contains the corresponding branch length. Pixels which are at
%   the root of their respective tree have MAP(x) = x and GAPS(x) =
%   inf.
%
%   [MAP,GAPS,E] = VL_QUICKSHIFT(I, KERNELSIZE, MAXDIST) also returns the estimate
%   of the density E.
%
%   [MAP,GAPS] = VL_QUICKSHIFT(I, KERNELSIZE) uses a default MAXDIST of 3 * KERNELSIZE.
%
%   Notes::
%     The distance between pixels is always measured in image
%     coordinates (not normalized), so the importance of the color
%     component should be weighted accordingly before calling this
%     function.
%
%   Options:
%
%   Verbose::
%     Toggles verbose output.
%
%   Medoid::
%     Run medoid shift instead of quick shift.
%
%   See also: VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function collocation_test ( )

%*****************************************************************************80
%
%% COLLOCATION_TEST tests the COLLOCATION library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 March 2015
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'COLLOCATION_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the COLLOCATION library.\n' );

  collocation_test01 ( );
  collocation_test02 ( );
  collocation_test03 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'COLLOCATION_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function I = not(I)

I.Z = ~I.Z;
I.name = ['not(' I.name ')'];

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
function [ objectIndex ] = GetScenario_randomUniform(varargin)

fprintf('[SCENARIO]\tGetting a normally distributed object scenario.\n');

% DEFAULT AGENT/OBSTACLE CONFIGURATION IN THIS EXAMPLE
defaultConfig = struct(...
    'file','scenario.mat',...
    'objects',[],...
    'waypointRadius',0.5,...
    'velocity',18,...
    'positionGain',10,...
    'velocityGain',1,...
    'poseGain',pi,...
    'noiseFactor',0,...
    'is3D',true,...
    'plot',false);

% Instanciate the scenario builder
SBinstance = scenarioBuilder();
% Parse user inputs
[inputConfig] = SBinstance.configurationParser(defaultConfig,varargin);

% Get the random object configuration
objectConfig = SBinstance.randomUniform(...
    'objects',numel(inputConfig.objects),...
    'positionGain',inputConfig.positionGain,...
    'velocityGain',inputConfig.velocityGain,...
    'poseGain',inputConfig.poseGain);

for index = 1:numel(inputConfig.objects)
    % Assign object
    objectIndex{index} = inputConfig.objects{index};
    % APPLY GLOBAL STATE VARIABLES
    if inputConfig.is3D
        objectIndex{index}.SetGLOBAL('position',objectConfig.positions(:,index) + inputConfig.noiseFactor*randn(3,1);
        objectIndex{index}.SetGLOBAL('velocity',objectConfig.velocities(:,index) + inputConfig.noiseFactor*randn(3,1);
        objectIndex{index}.SetGLOBAL('quaternion',objectConfig.quaternions(:,index));
    else
        objectIndex{index}.SetGLOBAL('position',[objectConfig.positions(1:2,index);0] + inputConfig.noiseFactor*[randn(2,1);0];
        objectIndex{index}.SetGLOBAL('velocity',[objectConfig.velocities(1:2,index);0] + inputConfig.noiseFactor*[randn(2,1);0];
        eta = OMAS_geometry.quaternionToEulers(objectConfig.quaternions(:,index));
        objectIndex{index}.SetGLOBAL('quaternion',OMAS_geometry.eulersToQuaternion([eta(3);0;0]));
    end
end
% Get the random object configuration
waypointConfig = SBinstance.randomUniform(...
    'objects',numel(inputConfig.objects),...
    'positionGain',inputConfig.positionGain,...
    'velocityGain',inputConfig.velocityGain,...
    'poseGain',inputConfig.poseGain);

fprintf('[SCENARIO]\tAssigning waypoints to agents.\n');
waypointSet = cell(size(inputConfig.objects));
for index = 1:numel(inputConfig.objects)
    % WAYPOINT SET A
    waypointSet{index} = waypoint('radius',inputConfig.waypointRadius,'priority',1,'name',sprintf('WP-%s',objectIndex{index}.name));
    % APPLY GLOBAL STATE VARIABLES
    if inputConfig.is3D
        waypointSet{index}.SetGLOBAL('position',waypointConfig.positions(:,index) + inputConfig.noiseFactor*randn(3,1));
        waypointSet{index}.SetGLOBAL('velocity',waypointConfig.velocities(:,index));
        waypointSet{index}.SetGLOBAL('quaternion',waypointConfig.quaternions(:,index));
    else
        waypointSet{index}.SetGLOBAL('position',[waypointConfig.positions(1:2,index);0] + inputConfig.noiseFactor*[randn(2,1);0]);
        waypointSet{index}.SetGLOBAL('velocity',[waypointConfig.velocities(1:2,index);0] + inputConfig.noiseFactor*[randn(2,1);0]);
        eta = OMAS_geometry.quaternionToEulers(waypointConfig.quaternions(:,index));
        waypointSet{index}.SetGLOBAL('quaternion',OMAS_geometry.eulersToQuaternion([eta(3);0;0]));
    end
    waypointSet{index} = waypointSet{index}.CreateAgentAssociation(objectIndex{index});  % Create waypoint with association to agent
end
% Concatinate the objects and way-points
objectIndex = [objectIndex;waypointSet];

% PLOT THE SCENE
if inputConfig.plot
    SBinstance.plotObjectIndex(objectIndex);
end
% CLEAR THE REMAINING VARIABLES
clearvars -except objectIndex
end

repo: OpenMAS


text:
function [model, grChek] = vargplvmOptimise(model, display, iters, varargin)

% VARGPLVMOPTIMISE Optimise the VARGPLVM.
% FORMAT
% DESC takes a given GP-LVM model structure and optimises with
% respect to parameters and latent positions.
% ARG model : the model to be optimised.
% ARG display : flag dictating whether or not to display
% optimisation progress (set to greater than zero) (default value 1).
% ARG iters : number of iterations to run the optimiser
% for (default value 2000).
% RETURN model : the optimised model.
%
% SEEALSO : vargplvmCreate, vargplvmLogLikelihood,
% vargplvmLogLikeGradients, vargplvmObjective, vargplvmGradient
%
% COPYRIGHT : Michalis K. Titsias, 2009
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006

% VARGPLVM

if nargin < 3
    iters = 2000;
    if nargin < 2
        display = 1;
    end
end

options = optOptions;
params = modelExtractParam(model);
if length(varargin) == 2
    if strcmp(varargin{1}, 'gradcheck')
        assert(islogical(varargin{2}));
        %options(9) = varargin{2};
        doGradchek = varargin{2};
        if doGradchek
            [gradient, delta] = feval('gradchek', params, @vargplvmObjective, @vargplvmGradient, model);
            deltaf = gradient - delta;
            d=norm(deltaf - gradient)/norm(gradient + deltaf);
            d1=norm(deltaf - gradient,1)/norm(gradient + deltaf,1);
            grRatio = sum(abs(gradient ./ deltaf)) / length(deltaf);
            fprintf(1,' Norm1 difference: %d\n Norm2 difference: %d\n Ratio: %d\n',d1,d, grRatio);
            grChek = {delta, gradient, deltaf, d, d1};
        else
            grChek = [];
        end
    end
end

options(2) = 0.1*options(2);
options(3) = 0.1*options(3);

if display
    options(1) = 1;
 %   if length(params) <= 100
 %       options(9) = 1;
 %   end
end
options(14) = iters;

if iters > 0
    if isfield(model, 'optimiser') && ~isa(model.optimiser, 'function_handle')
        if isfield(model, 'optimiser')
            optim = str2func(model.optimiser);
        else
            optim = str2func('scg');
        end

        if strcmp(func2str(optim), 'optimiMinimize')
            % Carl Rasmussen's minimize function
            params = optim('vargplvmObjectiveGradient', params, options, model);
        elseif strcmp(func2str(optim), 'scg2')
            % NETLAB style optimization with a slight modification so that an
            % objectiveGradient can be used where applicable, in order to re-use
            % precomputed quantities.
            params = optim('vargplvmObjectiveGradient', params,  options,  'vargplvmGradient', model);
        else
            % NETLAB style optimization.
            params = optim('vargplvmObjective', params,  options,  'vargplvmGradient', model);
        end
    elseif isfield(model, 'optimiser') && isa(model.optimiser, 'function_handle')
        f = fcnchk(model.optimiser);
        params = f(model);
    else
        error('vargplvmOptimise: Invalid optimiser setting.');
    end
    %model = vargplvmExpandParam(model, params);
    model = modelExpandParam(model, params);

    % Check if SNR of the optimised model is reasonable (ortherwise a
    % bad local minimum might have been found)
    if isfield(model, 'throwSNRError')
        svargplvmCheckSNR({vargplvmShowSNR(model)}, [], [], model.throwSNRError);
    else
        svargplvmCheckSNR({vargplvmShowSNR(model)});
    end
end

author: SheffieldML

repo: vargplvm


text:
function [m] = dm2m(dm)
% Convert length from decimeters to meters.
% Chad A. Greene 2012
m = dm/10;

author: Sable

repo: mcbench-benchmarks


text:
%function feat = rgb_fun(filename, fun)
% apply function handler fun separately to all rgb channel of image located
% in filename
function feat = yuv_fun(filename, fun)
rgbx = imread(filename);

feat = [];
if (numel(size(rgbx)) ~= 3)
    rgbx = repmat(rgbx,[1,1,3]);
end
yuvx = rgb2yuv(rgbx);
feat = [];
for c = 1:3
    cx = single(squeeze(yuvx(:,:,c)));
    feat = [feat, fun(cx)];
end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function h = matrix_decomp

%% Problem instance settings

m = 20;
n = 50;
use_cvx = 1; % set to 0 for larger instances

%% Problem data

s = RandStream.create('mt19937ar','seed',5489);
RandStream.setDefaultStream(s);

N = 3;
r = 4;

L = randn(m,r) * randn(r,n);    % low rank
S = sprandn(m,n,0.05);          % sparse
S(S ~= 0) = 20*binornd(1,0.5,nnz(S),1)-10;
V = 0.01*randn(m,n);            % noise

A = S + L + V;

g2_max = norm(A(:),inf);
g3_max = norm(A);
g2 = 0.15*g2_max;
g3 = 0.15*g3_max;

%% CVX

if use_cvx
    tic;

    cvx_begin
        cvx_precision low
        variables X_1(m,n) X_2(m,n) X_3(m,n)
        minimize(0.5*square_pos(norm(X_1,'fro')) + g2*norm(X_2(:),1) + g3*norm_nuc(X_3))
        subject to
            X_1 + X_2 + X_3 == A;
    cvx_end

    h.cvx_toc = toc;
    h.p_cvx = cvx_optval;
    h.X1_cvx = X_1;
    h.X2_cvx = X_2;
    h.X3_cvx = X_3;

    X_2(abs(X_2) < 1e-4) = 0;
    rhat = sum(svd(X_3) > 1e-4);
    fprintf('CVX (vs true):\n');
    fprintf('|V| = %.2f;  |X_1| = %.2f\n', norm(V, 'fro'), norm(X_1,'fro'));
    fprintf('nnz(S) = %d; nnz(X_2) = %d\n', nnz(S), nnz(X_2));
    fprintf('rank(L) = %d; rank(X_3) = %d\n', rank(L), rhat);
end

%% ADMM

MAX_ITER = 100;
ABSTOL   = 1e-4;
RELTOL   = 1e-2;

tic;

lambda = 1;
rho = 1/lambda;

X_1 = zeros(m,n);
X_2 = zeros(m,n);
X_3 = zeros(m,n);
z   = zeros(m,N*n);
U   = zeros(m,n);

fprintf('\n%3s\t%10s\t%10s\t%10s\t%10s\t%10s\n', 'iter', ...
    'r norm', 'eps pri', 's norm', 'eps dual', 'objective');

for k = 1:MAX_ITER

    B = avg(X_1, X_2, X_3) - A./N + U;

    % x-update
    X_1 = (1/(1+lambda))*(X_1 - B);
    X_2 = prox_l1(X_2 - B, lambda*g2);
    X_3 = prox_matrix(X_3 - B, lambda*g3, @prox_l1);

    % (for termination checks only)
    x = [X_1 X_2 X_3];
    zold = z;
    z = x + repmat(-avg(X_1, X_2, X_3) + A./N, 1, N);

    % u-update
    U = B;

    % diagnostics, reporting, termination checks
    h.objval(k)   = objective(X_1, g2, X_2, g3, X_3);
    h.r_norm(k)   = norm(x - z,'fro');
    h.s_norm(k)   = norm(-rho*(z - zold),'fro');
    h.eps_pri(k)  = sqrt(m*n*N)*ABSTOL + RELTOL*max(norm(x,'fro'), norm(-z,'fro'));
    h.eps_dual(k) = sqrt(m*n*N)*ABSTOL + RELTOL*sqrt(N)*norm(rho*U,'fro');

    if k == 1 || mod(k,10) == 0
        fprintf('%4d\t%10.4f\t%10.4f\t%10.4f\t%10.4f\t%10.2f\n', k, ...
            h.r_norm(k), h.eps_pri(k), h.s_norm(k), h.eps_dual(k), h.objval(k));
    end

    if h.r_norm(k) < h.eps_pri(k) && h.s_norm(k) < h.eps_dual(k)
         break;
    end

end

h.admm_toc = toc;
h.admm_iter = k;
h.X1_admm = X_1;
h.X2_admm = X_2;
h.X3_admm = X_3;

fprintf('\nADMM (vs true):\n');
fprintf('|V| = %.2f;  |X_1| = %.2f\n', norm(V, 'fro'), norm(X_1,'fro'));
fprintf('nnz(S) = %d; nnz(X_2) = %d\n', nnz(S), nnz(X_2));
fprintf('rank(L) = %d; rank(X_3) = %d\n', rank(L), rank(X_3));

if use_cvx
    fprintf('\nADMM vs CVX solutions (in Frobenius norm):\n');
    fprintf('X_1: %.2e; X_2: %.2e; X_3: %.2e\n', ...
        norm(h.X1_cvx - X_1,'fro'), norm(h.X2_cvx - X_2,'fro'), norm(h.X3_cvx - X_3,'fro'));
end

end

function x = avg(varargin)
    N = length(varargin);
    x = 0;
    for k = 1:N
        x = x + varargin{k};
    end
    x = x/N;
end

function p = objective(X_1, g_2, X_2, g_3, X_3)
    p = norm(X_1,'fro').^2 + g_2*norm(X_2(:),1) + g_3*norm(svd(X_3),1);
end

author: andrewssobral

repo: lrslibrary


text:
function timer_tictoc_test01 ( )

%*****************************************************************************80
%
%% TEST01 times the RAND routine.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 June 2006
%
%  Author:
%
%    John Burkardt
%
  n_log_min = 10;
  n_log_max = 20;
  n_min = 2^n_log_min;
  n_max = 2^n_log_max;
  n_rep = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01\n' );
  fprintf ( 1, '  Time the MATLAB RAND routine:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    x = rand(n,1);\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Data vectors will be of minimum size %d\n', n_min );
  fprintf ( 1, '  Data vectors will be of maximum size %d\n', n_max );
  fprintf ( 1, '  Number of repetitions of the operation: %d\n', n_rep );

  for i_rep = 1 : n_rep

    for n_log = n_log_min : n_log_max

      n = 2^n_log;

      tic;

      x = rand(n,1);

      delta(n_log,i_rep) = toc;

    end

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST01 Results:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Vector Size  Rep #1        Rep #2        Rep #3        ' );
  fprintf ( 1, 'Rep #4        Rep #5\n' );
  fprintf ( 1, '\n' );

  for n_log = n_log_min : n_log_max
    n = 2^n_log;
    fprintf ( 1, '%10d', n );
    for j = 1 : n_rep
      fprintf ( '%14f', delta(n_log,j) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


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

%%
clear; close all; clc;

%% Syntax
% |[G]=patchFeatureDetect(F,V,a);|

%% Description
% This function detects surface features or groups for the input patch data
% defined by the faces F and the vertices V. Patch elements for which the
% dihedral angle is lower than a are grouped together.
% This function is useful for detecting sets of faces from for instance
% importance CAD geometry, e.g. top faces, side faces etc can be
% automatically detected provided that their boundaries have a significant
% enough angle change.

%%
% Plot settings
fontSize=15;

%%
% Create example data

testCase=4;
switch testCase
    case 1
        boxDim=[4 5 6]; %Width in each direction
        pointSpacing=1; %Desired point spacing
        [F,V]=triBox(boxDim,pointSpacing);
    case 2
        boxDim=[4 5 6]; %Width in each direction
        boxEl=[4 5 6]; %Desired number of elements
        [F,V]=quadBox(boxDim,boxEl);
    case 3
        defaultFolder = fileparts(fileparts(mfilename('fullpath')));
        pathName=fullfile(defaultFolder,'data','libSurf');
        dataStruct=load(fullfile(pathName,'enginePart_p1.mat'));
        F=dataStruct.F;
        V=dataStruct.V;
    case 4
        defaultFolder = fileparts(fileparts(mfilename('fullpath')));
        pathName=fullfile(defaultFolder,'data','libSurf');
        dataStruct=load(fullfile(pathName,'sprocket.mat'));
        F=dataStruct.F;
        V=dataStruct.V;
end

%%
% Use |patchFeatureDetect| to detect features in patch

%Angular threshold in radians
a=(45/180)*pi;

%Detect surface features
G=patchFeatureDetect(F,V,a);

%%

% Plotting model
cFigure;
subplot(1,2,1); hold on;
title('Input surface','FontSize',fontSize);
gpatch(F,V,'w','k',1);
axisGeom(gca,fontSize); camlight headlight;

subplot(1,2,2); hold on;
title('Surface feature labels','FontSize',fontSize);
gpatch(F,V,G,'k',1);
axisGeom(gca,fontSize); camlight headlight;
colormap gjet; icolorbar;
gdrawnow;

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

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [v,d,z]=soundspeed(t,p,m,g)
%SOUNDSPEED gives the speed of sound, density of air and acostuc impedance as a function of temp & pressure [V,D,Z]=(T,P,M,G)
%
%  Inputs:  T        air temperature in Celsius  [20 deg C]
%           P        air pressure [1 atm]
%           M        average molecular weight of air [0.0289644 kg/mol]
%           G        adiabatic constant for air [1.4]
%
% Outputs:  V        is the speed of sound in m/s
%           D        density of air in kg/m^3
%           Z        characteristic impedance of air Pa.s/m

% Notes: (1) Sound pressure is often measured in dB (SPL) relative to 20uPa [20*log10(p/p0)]
%            Sound pressure is inversely proportional to distance.
%        (2) Sound intensity is often measured indB relative to pW/m^2,[10*log10(J*10^12)]
%            Intensity is inversely proportional to distance squared.
%        (3) Intensity * impedance = pressure^2, so with the default values, 1 pW/m^2 = 20.33 uPa
%            So: X dB (SPL) = X-93.98 dB (Pa) = X-0.14 dB (pW/m^2) =  X-120.14 dB (W/m^2)
%        (4) The default air pressure (which does not affect sound speed) in various units is:
%            1 atm = 101325 Pa = 1.01325 bar = 1.0332 at = 760 torr = 14.696 psi

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

if nargin<4
    g=1.4;
    if nargin<3
        m=0.0289644;      % gm/mol
        if nargin<2
            p=1;
            if nargin<1
                t=20;
            end
        end
    end
end
p=p*101325; % convert pressure: atm to pascal
k=t+273.15;  % absolute temperature
r=8.3144;  % J/(mol K) universal gas constant
d=p*m/(r*k);
v=sqrt(g*r*k/m);
z=v*d;

author: decouples


text:
function pass = test_chebcoeffs3(pref)
% Test chebcoeffs3 and coeffs3.

if ( nargin == 0)
    pref = chebfunpref;
end
tol = 1000*pref.cheb3Prefs.chebfun3eps;

% Create a rank-(3,3,3) function:
m = 8;
n = 10;
p = 5;
Tm = chebpoly(m);
Tn = chebpoly(n);
Tp = chebpoly(p);
f = outerProd(Tm, Tm, Tm) + outerProd(Tn, Tn, Tn) + outerProd(Tp, Tp, Tp);

% Create expected tensor of coefficients for f:
fCoeffs = chebcoeffs3(f);
coeffsExact = zeros(m+1, n+1, p+1);
coeffsExact(m+1, m+1, m+1) = 1;
coeffsExact(n+1, n+1, n+1) = 1;
coeffsExact(p+1, p+1, p+1) = 1;
pass(1) = norm(fCoeffs(:) - coeffsExact(:)) < tol;

% Check the same with coeffs3 instead of chebcoeffs3:
fCoeffs = coeffs3(f);
coeffsExact = zeros(m+1, n+1, p+1);
coeffsExact(m+1, m+1, m+1) = 1;
coeffsExact(n+1, n+1, n+1) = 1;
coeffsExact(p+1, p+1, p+1) = 1;
pass(2) = norm(fCoeffs(:) - coeffsExact(:)) < tol;

% check the reverse:
C = chebfun3.vals2coeffs(chebpolyval3(f));
pass(3) = norm(C(:) - coeffsExact(:)) < tol;

% Check that it works for multiple outputs also:
[core, cols_coeffs, rows_coeffs, tubes_coeffs] = chebcoeffs3(f);
fCoeffs = chebfun3.txm(chebfun3.txm(chebfun3.txm(core, ...
        cols_coeffs, 1), rows_coeffs, 2), tubes_coeffs, 3);
coeffsExact = zeros(m+1, n+1, p+1);
coeffsExact(m+1, m+1, m+1) = 1;
coeffsExact(n+1, n+1, n+1) = 1;
coeffsExact(p+1, p+1, p+1) = 1;
pass(4) = norm(fCoeffs(:) - coeffsExact(:)) < tol;

% Check that coeffs3 also works for multiple outputs also:
[core, cols_coeffs, rows_coeffs, tubes_coeffs] = coeffs3(f);
fCoeffs = chebfun3.txm(chebfun3.txm(chebfun3.txm(core, ...
        cols_coeffs, 1), rows_coeffs, 2), tubes_coeffs, 3);
coeffsExact = zeros(m+1, n+1, p+1);
coeffsExact(m+1, m+1, m+1) = 1;
coeffsExact(n+1, n+1, n+1) = 1;
coeffsExact(p+1, p+1, p+1) = 1;
pass(5) = norm(fCoeffs(:) - coeffsExact(:)) < tol;

end

author: chebfun

repo: chebfun


text:
function codingLength = total_coding_length(data, labels, epsilon, affine)

% total_coding_length
%
%   compute the total number of bits needed to code a set of data vectors
%   given a grouping, and an error tolerance.
%
% Inputs:
%   data - a matrix where each column is a data vector.
%   labels - a list of integer labels that partitions the data into groups
%   epsilon - the error tolerance
%   affine - if true, compute using the affine coding length function.
%               otherwise use the linear coding length function
%
% Outputs:
%   codingLength - the total number of bits needed to code this set of data
%
% Dependencies:
%   coding_length
%
% Sep. '07  Shankar Rao -- srrao@uiuc.edu

% Copyright 2007, University of Illinois. All rights reserved.

codingLength = 0;
groupCount = max(labels);
sampleCount = length(labels);
% Make sure outlier indices are counted too
% modified 10-24-07
setOfIndices = unique(labels);
for groupIndex = setOfIndices(:)',
    codingLength = codingLength + coding_length(data(:, labels == groupIndex), epsilon^2, sampleCount, affine);
end

author: SuTanTank

repo: VideoStitchingViaShakinessRemoving


text:
function hb_to_msm_test ( )

%*****************************************************************************80
%
%% HB_TO_MSM_TEST tests the HB_TO_MSM library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 August 2008
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HB_TO_MSM_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the HB_TO_MSM library.\n' );

  hb_to_msm_test01 ( );
  hb_to_msm_test02 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HB_TO_MSM_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
function y = quantile1(x,p)
%QUANTILE1 Quantile of a vector.
%   Y = PRCTILE(X,P) returns percentiles of the values in X.  P is a scalar
%   or a vector of percent values.  When X is a vector, Y is the same size
%   as P, and Y(i) contains the P(i)-th percentile.  When X is a matrix,
%   the i-th row of Y contains the P(i)-th percentiles of each column of X.
%   For N-D arrays, PRCTILE operates along the first non-singleton
%   dimension.
%
%   Percentiles are specified using percentages, from 0 to 100.  For an N
%   element vector X, PRCTILE computes percentiles as follows:
%      1) The sorted values in X are taken as the 100*(0.5/N), 100*(1.5/N),
%         ..., 100*((N-0.5)/N) percentiles.
%      2) Linear interpolation is used to compute percentiles for percent
%         values between 100*(0.5/N) and 100*((N-0.5)/N)
%      3) The minimum or maximum values in X are assigned to percentiles
%         for percent values outside that range.
%
%   PRCTILE treats NaNs as missing values, and removes them.
%
%   Examples:
%      y = prctile(x,50); % the median of x
%      y = prctile(x,[2.5 25 50 75 97.5]); % a useful summary of x
%
%   See also IQR, MEDIAN, NANMEDIAN, QUANTILE.

%   Copyright 1993-2016 The MathWorks, Inc.

% If X is empty, return all NaNs.
if isempty(x)
    y = nan(size(p),'like',x);
else
    % Drop X's leading singleton dims, and combine its trailing dims.  This
    % leaves a matrix, and we can work along columns.
    x = x(:);

    x = sort(x,1);
    n = sum(~isnan(x), 1); % Number of non-NaN values

    if isequal(p,0.5) % make the median fast
        if rem(n,2) % n is odd
            y = x((n+1)/2,:);
        else        % n is even
            y = (x(n/2,:) + x(n/2+1,:))/2;
        end
    else
        r = p*n;
        k = floor(r+0.5); % K gives the index for the row just before r
        kp1 = k + 1;      % K+1 gives the index for the row just after r
        r = r - k;        % R is the ratio between the K and K+1 rows

        % Find indices that are out of the range 1 to n and cap them
        k(k<1 | isnan(k)) = 1;
        kp1 = bsxfun( @min, kp1, n );

        % Use simple linear interpolation for the valid percentages
        y = (0.5+r).*x(kp1,:)+(0.5-r).*x(k,:);

        % Make sure that values we hit exactly are copied rather than interpolated
        exact = (r==-0.5);
        if any(exact)
            y(exact,:) = x(k(exact),:);
        end

        % Make sure that identical values are copied rather than interpolated
        same = (x(k,:)==x(kp1,:));
        if any(same(:))
            x = x(k,:); % expand x
            y(same) = x(same);
        end

    end

end

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc