text:
classdef LaunchVehicleState < matlab.mixin.SetGet & matlab.mixin.Copyable
    %LaunchVehicleState Summary of this class goes here
    %   Detailed explanation goes here

    properties
        lv LaunchVehicle

        e2TConns EngineToTankConnState
        t2TConns TankToTankConnState
        holdDownEnabled(1,1) logical = false
    end

    properties(Transient)
        cachedEngines LaunchVehicleEngine
        cachedConnTanks cell
    end

    properties(Constant)
        emptyTankArr = LaunchVehicleTank.empty(1,0);
    end

    methods
        function obj = LaunchVehicleState(lv)
            obj.lv = lv;
            obj.clearCachedConnEnginesTanks();
        end

        %Engine to Tank Connections
        function addE2TConnState(obj, newConnState)
            obj.e2TConns(end+1) = newConnState;

            obj.clearCachedConnEnginesTanks();
        end

        function removeE2TConnStateForConn(obj, conn)
            ind = [];
            for(i=1:length(obj.e2TConns)) %#ok<*NO4LP>
                if(obj.e2TConns(i).conn == conn)
                    ind = i;
                    break;
                end
            end

            if(not(isempty(ind)))
                obj.e2TConns(ind) = [];
            end

            obj.clearCachedConnEnginesTanks();
        end

        function tanks = getTanksConnectedToEngine(obj, engine)
            bool = obj.cachedEngines == engine;
            if(any(bool))
                tanks = obj.cachedConnTanks{bool};
            else
                tanks = obj.emptyTankArr;

                connStates = obj.e2TConns([obj.e2TConns.active] == true);
                connections = [connStates.conn];

                if(not(isempty(connections)))
                    connections = connections([connections.engine] == engine);
                    tanks = [connections.tank];
                end

                if(isempty(tanks))
                    tanks = obj.emptyTankArr;
                end

                obj.cachedEngines(end+1) = engine;
                obj.cachedConnTanks(end+1) = {tanks};
            end
        end

        function tanks = getTanksConnectedToEngineForStage(obj, engine, stage)
%             tanks = obj.emptyTankArr;

            engineTanks = obj.getTanksConnectedToEngine(engine);
%             for(i=1:length(engineTanks))
%                 if(engineTanks(i).stage == stage)
%                     tanks(end+1) = engineTanks(i); %#ok<AGROW>
%                 end
%             end
            if(not(isempty(engineTanks)))
                bool = [engineTanks.stage] == stage;
                tanks = engineTanks(bool);
            else
                tanks = obj.emptyTankArr;
            end
        end

        %Tank To Tank Connections
        function addT2TConnState(obj, newConnState)
            obj.t2TConns(end+1) = newConnState;
        end

        function removeT2TConnStateForConn(obj, conn)
            ind = [];
            for(i=1:length(obj.t2TConns)) %#ok<*NO4LP>
                if(obj.t2TConns(i).conn == conn)
                    ind = i;
                    break;
                end
            end

            if(not(isempty(ind)))
                obj.t2TConns(ind) = [];
            end
        end

        function t2TConnState = getTank2TankConnStateForConn(obj, conn)
            ind = [];
            for(i=1:length(obj.t2TConns)) %#ok<*NO4LP>
                if(obj.t2TConns(i).conn == conn)
                    ind = i;
                    break;
                end
            end

            if(not(isempty(ind)))
                t2TConnState = obj.t2TConns(ind);
            end
        end

        function tanks = getTanksWithActiveTankToTankConnectionsForStage(obj, stageState)
            tanks = LaunchVehicleState.emptyTankArr;

            for(i=1:length(obj.t2TConns))
                connState = obj.t2TConns(i);

                if(connState.active)
                    conn = connState.conn;

                    if((not(isempty(conn.srcTank)) && conn.srcTank.stage == stageState.stage))
                        tanks(end+1) = conn.srcTank; %#ok<AGROW>
                    end

                    if((not(isempty(conn.tgtTank)) && conn.tgtTank.stage == stageState.stage))
                        tanks(end+1) = conn.tgtTank; %#ok<AGROW>
                    end
                end
            end

            if(length(tanks) > 1)
                tanks = unique(tanks);
            end
        end

        %Misc
        function newLvState = deepCopy(obj)
            newLvState = obj.copy();

            for(i=1:length(obj.e2TConns))
                newLvState.e2TConns(i) = obj.e2TConns(i).deepCopy();
            end

            for(i=1:length(obj.t2TConns))
                newLvState.t2TConns(i) = obj.t2TConns(i).deepCopy();
            end

            newLvState.clearCachedConnEnginesTanks();
        end

        function clearCachedConnEnginesTanks(obj)
            obj.cachedEngines = LaunchVehicleEngine.empty(1,0);
            obj.cachedConnTanks = {};
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function M = centeredmatrixfactory(m, n, rows_or_cols)
% Linear manifold struct. for optimization over matrices with centered cols
%
% function M = centeredmatrixfactory(m, n)
% function M = centeredmatrixfactory(m, n, 'cols')
% function M = centeredmatrixfactory(m, n, 'rows')
%
% Returns M, a structure for Manopt describing the Euclidean space of
% m-by-n matrices whose columns sum to zero (or whose rows sum to zero,
% if 'rows' is passed as last input).
%
% The metric is the standard Frobenius distance and associated trace inner
% product. Matrices on M, denoted by X, have size mxn and obey
% X*ones(n, 1) = 0 (centered columns) or ones(1, m)*X = 0 (centered rows).
%
% See also: euclideanfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, July 3, 2015.
% Contributors:
% Change log:

    if ~exist('rows_or_cols', 'var') || isempty(rows_or_cols)
        rows_or_cols = 'cols';
    end

    % Define a centering operator: it subtracts the mean column or row.
    switch lower(rows_or_cols)
        case 'cols'
            center = @(X) bsxfun(@minus, X, mean(X, 2));
            M.dim = @() m*n - m;
        case 'rows'
            center = @(X) bsxfun(@minus, X, mean(X, 1));
            M.dim = @() m*n - n;
        otherwise
            error('The third input must be either ''rows'' or ''cols''.');
    end

    % This is a non-standard function to have in a Manopt manifold.
    % It is included because it might be helpful in some situations.
    M.center = center;

    M.name = @() sprintf('Space of size %d x %d matrices with centered %s', ...
                         m, n, lower(rows_or_cols));

    M.inner = @(x, d1, d2) d1(:).'*d2(:);

    M.norm = @(x, d) norm(d, 'fro');

    M.dist = @(x, y) norm(x-y, 'fro');

    M.typicaldist = @() sqrt(M.dim());

    M.proj = @(X, U) center(U);

    M.egrad2rgrad = M.proj;

    M.ehess2rhess = @(x, eg, eh, d) center(eh);

    M.tangent = @(x, d) d;

    M.exp = @exp;
    function y = exp(x, d, t)
        if nargin == 3
            y = x + t*d;
        else
            y = x + d;
        end
    end

    M.retr = M.exp;

	M.log = @(x, y) y-x;

    M.hash = @(x) ['z' hashmd5(x(:))];

    M.randvec = @(X) randvec();
    function U = randvec()
        U = center(randn(m, n));
        U = U / norm(U, 'fro');
    end

    M.rand = @() center(randn(m, n));

    M.lincomb = @matrixlincomb;

    M.zerovec = @(x) zeros(m, n);

    M.transp = @(x1, x2, d) d;

    M.pairmean = @(x1, x2) .5*(x1+x2);

    M.vec = @(x, u_mat) u_mat(:);
    M.mat = @(x, u_vec) reshape(u_vec, [m, n]);
    M.vecmatareisometries = @() true;

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function len = polygonLength(poly, varargin)
%POLYGONLENGTH Perimeter of a polygon.
%
%   L = polygonLength(POLYGON);
%   Computes the boundary length of a polygon. POLYGON is given by a N-by-2
%   array of vertices.
%
%   Example
%     % Perimeter of a circle approximation
%     poly = circleToPolygon([0 0 1], 200);
%     polygonLength(poly)
%     ans =
%         6.2829
%
%   See also
%   polygons2d, polygonCentroid, polygonArea, drawPolygon, polylineLength
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2005-05-11
% Copyright 2005-2022 INRA - TPV URPOI - BIA IMASTE

% If first argument is a cell array, this is a multi-polygon, and we simply
% add the lengths of individual polygons
if iscell(poly)
    len = 0;
    for i = 1:length(poly)
        len = len + polygonLength(poly{i});
    end
    return;
end

% case of a polygon given as two coordinate arrays
if nargin == 2
    poly = [poly varargin{1}];
end

% check there are enough points
if size(poly, 1) < 2
    len = 0;
    return;
end

% compute length
if size(poly, 2) == 2
    % polygon in dimension 2 (classical case)
    dp = diff(poly([1:end 1], :), 1, 1);
    len = sum(hypot(dp(:, 1), dp(:, 2)));
else
    % polygon of larger dimension
    len = sum(sqrt(sum(diff(poly([2:end 1], :), 1, 1).^2, 2)));
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function transition_net_multi(varargin)

opts.expDir = fullfile('../data','DCFNet-net-7') ;
opts = vl_argparse(opts, varargin);

modelPath = dir(fullfile(opts.expDir, 'net-epoch-*0.mat'));
modelPath = sort({modelPath.name});

[~,out_name,~] = fileparts(opts.expDir);

for i = 1:numel(modelPath)
    load(fullfile(opts.expDir, modelPath{i}));
    net = dagnn.DagNN.loadobj(net) ;
    net = deployDCFNet(net);
    if ~exist('../model_multi', 'dir')
        mkdir('../model_multi');
    end
    save(fullfile('../model_multi', [out_name num2str(i, 'epoch_%d.mat')]), 'net') ;
end
end

repo: UDT


text:
% Image Reconstruction With Directional Sensors Example
%
% This example demonstrates how the directionality of sensor elements can
% give rise to artefacts in time reversal photoacoustic image
% reconstruction. It builds on the Sensor Element Directivity in 2D and 2D
% Time Reversal Reconstruction For A Line Sensor examples.
%
% author: Bradley Treeby & Ben Cox
% date: 18th January 2010
% last update: 24th August 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

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

clear all

% =========================================================================
% SIMULATION
% =========================================================================

% create the computational grid
PML_size = 20;          % size of the PML in grid points
Nx = 128 - 2*PML_size;  % number of grid points in the x (row) direction
Ny = 256 - 2*PML_size;  % number of grid points in the y (column) direction
dx = 0.1e-3;            % grid point spacing in the x direction [m]
dy = 0.1e-3;            % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy);

% define the properties of the propagation medium
medium.sound_speed = 1500;	% [m/s]

% create initial pressure distribution using makeDisc
disc_magnitude = 5; % [au]
disc_x_pos = 60;    % [grid points]
disc_y_pos = 140;  	% [grid points]
disc_radius = 5;    % [grid points]
disc_2 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

disc_x_pos = 30;    % [grid points]
disc_y_pos = 110; 	% [grid points]
disc_radius = 8;    % [grid points]
disc_1 = disc_magnitude*makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

% smooth the initial pressure distribution and restore the magnitude
p0 = smooth(kgrid, disc_1 + disc_2, true);

% assign to the source structure
source.p0 = p0;

% define a four-sided, square sensor
sensor.mask = zeros(kgrid.Nx, kgrid.Ny);
sensor.mask(1, :) = 1;
sensor.mask(end, :) = 1;
sensor.mask(:, 1) = 1;
sensor.mask(:, end) = 1;

% create the time array
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);

% set the input arguements
input_args = {'PMLInside', false, 'PMLSize', PML_size, 'PlotPML', false, 'Smooth', false};

% run the simulation for omnidirectional detector elements
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

% define the directionality of the sensor elements
sensor.directivity_angle = zeros(kgrid.Nx, kgrid.Ny);
sensor.directivity_angle(1, :) = 0;    	 % max sensitivity in x direction
sensor.directivity_angle(end, :) = 0;  	 % max sensitivity in x direction
sensor.directivity_angle(:, 1) = pi/2;   % max sensitivity in y direction
sensor.directivity_angle(:, end) = pi/2; % max sensitivity in y direction

% define the directivity size
sensor.directivity_size = 20*kgrid.dx;

% run the simulation with directional elements
sensor_data_directional = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

% reset the initial pressure
source.p0 = 0;

% assign the time reversal data for the omnidirectional case
sensor.time_reversal_boundary_data = sensor_data;

% run the time reversal reconstruction
p0_recon = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

% assign the time reversal data for the directional case
sensor.time_reversal_boundary_data = sensor_data_directional;

% run the time reversal reconstruction with directional elements
p0_recon_directional = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

% =========================================================================
% VISUALISATION
% =========================================================================

% plot the initial pressure and sensor distribution
figure;
imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, p0 + sensor.mask*disc_magnitude, [-disc_magnitude disc_magnitude]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
colorbar;

% plot the reconstructed initial pressure
figure;
imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, p0_recon, [-disc_magnitude disc_magnitude]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
colorbar;

% plot the reconstructed initial pressure with directivity
figure;
imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, p0_recon_directional, [-disc_magnitude disc_magnitude]);
colormap(getColorMap);
ylabel('x-position [mm]');
xlabel('y-position [mm]');
axis image;
colorbar;

% plot a profile for comparison
figure;
plot(kgrid.y_vec*1e3, p0(disc_x_pos, :), 'k-',...
    kgrid.y_vec*1e3, p0_recon(disc_x_pos, :), 'r--',...
    kgrid.y_vec*1e3, p0_recon_directional(disc_x_pos, :), 'b:');
xlabel('y-position [mm]');
ylabel('Pressure');
legend('True', 'Omnidirectional','Directional');
axis tight;
set(gca, 'YLim', [0 5.1]);

repo: Aero-matlab


text:
function install_paths()
%THE SCRIPT to install the sltoolbox paths to matlab system

fp = mfilename('fullpath');
rootdir = fileparts(fp);

subfolders = { ...
    'ann'; ...
    'cluster'; ...
    'core'; ...
    'discrete'; ...
    'ExpDL'; ...
    'fileio'; ...
    'graph'; ...
    'imgproc'; ...
    'interp'; ...
    'kernel'; ...
    'learn'; ...
    'manifold'; ...
    'perfeval'; ...
    'regression'; ...
    'smallmat'; ...
    'stat'; ...
    'subspace'; ...
    'subspace_ex'; ...
    'tensor'; ...
    'text'; ...
    'utils'; ...
    'utils_ex'; ...
    'visualize'; ...
    'xmlkits'};

n = length(subfolders);
folderpaths = cell(1, n);

for i = 1 : n
    folderpaths{i} = fullfile(rootdir, subfolders{i});
    fprintf('Add path: %s\n', folderpaths{i});
end

addpath(folderpaths{:});
savepath;

disp('All paths have been added.');
disp(' ');  % a blank line

author: lmthang


text:
function Kc = kernelCenter(K)

% KERNELCENTER Attempts to Center Kernel Matrix
% FORMAT
% DESC returns a centered kernel matrix
% ARG kernel matrix
% RETURN Kc : The centered kernel
%
% SEEALSO :
%
% COPYRIGHT : Carl Henrik Ek, 2008

% KERN

if(nargin<1)
  error('To Few Arguments');
end
if(size(K,1)~=size(K,2))
  error('Kernel Not Square');
end

Kc = (eye(size(K)) - 1/size(K,1)*ones(size(K)))*K*(eye(size(K))-1/size(K,1)*ones(size(K)));

return;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [W,W0] = setWeight(W,p)
% Set the weight vectors

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

% This function is written by Tomoaki Takagi

    if 0 <= p && p <= 1
        %% Static approach
        W0 = [];
        M = size(W,2);
        % Distribution control of weight vector set
        TF = W < 1.0 / M;
        W(TF) = W(TF) * p * M;
        W(~TF) = 1.0 - (1.0 - W(~TF)) * (1.0 - p) * M / (M - 1);
    else
        %% Dynamic approach
        W0 = W;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function geometry_test2071 ( )

%*****************************************************************************80
%
%% TEST2071 tests TRIANGLE_POINT_DIST_2D;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
  ntest = 7;

  ptest = [ ...
     0.25,   0.25; ...
     0.75,   0.25; ...
     1.00,   1.00; ...
    11.00,   0.50; ...
     0.00,   1.00; ...
     0.50, -10.00; ...
     0.60,   0.60 ]';
  t = [ ...
    0.0, 1.0; ...
    0.0, 0.0; ...
    1.0, 0.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST2071\n' );
  fprintf ( 1, '  For a triangle in 2D,\n' );
  fprintf ( 1, '  TRIANGLE_POINT_DIST_2D computes the distance\n' );
  fprintf ( 1, '  to a point;\n' );

  r8mat_transpose_print ( 2, 3, t, '  Triangle vertices:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '       P       DIST\n' );
  fprintf ( 1, '\n' );

  for i = 1 : ntest

    p(1:2,1) = ptest(1:2,i);

    dist = triangle_point_dist_2d ( t, p );

    fprintf ( 1, '  %10f  %10f  %10f\n', p(1:2,1), dist );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%Sigmoid function
function[y] = logistic(x)
y = 1./(1 + exp(-x));

author: marianux

repo: ecg-kit

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


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%
% 	Problem 3- closed-loop discrete time Transfer Function

K=2;
Ts=-1

G=tf(0.1, [1 -0.5], Ts);

H = tf(0.5, [1 -0.1], Ts);

F = feedback(K*G, H)

author: Sable

repo: mcbench-benchmarks


text:
function [f_alias, cond_N] = sphArrayAliasLim(R, Nmic, maxN, mic_dirs_rad, mic_weights)
%SPHARRAYALIASLIM Get estimates of the aliasing limit of a spherical array
%
%   First estimate takes into account only the radius and a nominal order
%   that the array is expected to support, it is the simplest one and it
%   expresses the kR = maxN rule.
%   The second estimate is based on the number of microphones, and it can
%   be more relaxed than the first, if the nominal supported order is less
%   than maxN<floor(sqrt(Nmic)-1).
%   The third takes into account microphone numbers and directions, and it
%   is based on the orthogonality of the SH matrix for the microphone
%   positions expressed through the condition number.
%
%   Inputs:
%       R:      radius of spherical array
%       Nmic:   number of microphones
%       maxN:   expected maximum supported order of the array
%       mic_dirs_rad:   nMicx2 matrix of [azi elev] directions of the
%           microphones on the array
%       mic_weights: (optional) nMicx1 vector of weights used to improve
%           orthogonality of the SH transform
%
%   Outputs:
%       f_alias:    the alising limit estimates
%       cond_N:     the condition number of all orders up to one above
%           floor(sqrt(Nmic)-1)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SPHARRAYALIASLIM.M - 11/7/2013
% Archontis Politis, archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c = 343;

% conventional kR=N assumption
f_alias(1) = c*maxN/(2*pi*R);

% based on the floor of the number of microphones, uniform arrangement
f_alias(2) = c*floor(sqrt(Nmic)-1)/(2*pi*R);

% based on condition number of the SHT matrix
aziElev2aziIncl = @(dirs) [dirs(:,1) pi/2-dirs(:,2)];

maxOrder = ceil(sqrt(Nmic)-1);
if (nargin<5), mic_weights = []; end

cond_N = checkCondNumberSHT(maxOrder, aziElev2aziIncl(mic_dirs_rad), 'real', mic_weights);
% plot(cond_N)
trueMaxOrder = find(cond_N<10^4, 1, 'last')-1;
f_alias(3) = c*trueMaxOrder/(2*pi*R);

end

author: polarch

repo: Spherical-Array-Processing


text:
% Script to reproduce the experiments leading to the results provided in the
% Table 2 of the paper "Deep Scattering Spectrum" by J. Andén and S. Mallat.

% M=1 scattering

run_name = 'DSS_Table2_GTZAN_m1';

N=5*2^17;

src=gtzan_src('/path/to/gtzan');

fparam.filter_type = {'gabor_1d','morlet_1d'};
fparam.Q = [8 2];
fparam.J = T_to_J(8192,fparam);

options.M = 1;

Wop = wavelet_factory_1d(N, fparam, options);

feature_fun = {@(x)(format_scat(log_scat(renorm_scat(scat(x,Wop)))))};

db = prepare_database(src,feature_fun);
db.features = single(db.features);
db = svm_calc_kernel(db,'gaussian','square',1:2:size(db.features,2));

rs = RandStream.create('mt19937ar','Seed',floor(pi*1e9));
RandStream.setGlobalStream(rs);
[train_set{1}, test_set{1}] = create_partition([src.objects.class], 0.9);
for k = 2:10
	[train_set{k}, test_set{k}] = ...
		next_fold([src.objects.class], train_set{k-1}, test_set{k-1});
end

optt.kernel_type = 'gaussian';
optt.C = 2.^[0:4:8];
optt.gamma = 2.^[-16:4:-8];
optt.search_depth = 3;
optt.full_test_kernel = 1;

for k = 1:10
	[dev_err_grid,C_grid,gamma_grid] = ...
		svm_adaptive_param_search(db,train_set{k},[],optt);

	[dev_err(k),ind] = min(mean(dev_err_grid{end},2));
	C(k) = C_grid{end}(ind);
	gamma(k) = gamma_grid{end}(ind);

	optt1 = optt;
	optt1.C = C(k);
	optt1.gamma = gamma(k);

	model = svm_train(db,train_set{k},optt1);
	labels(:,k) = svm_test(db,model,test_set{k});
	err(k) = classif_err(labels(:,k),test_set{k},db.src);

	fprintf('dev err = %f, test err = %f\n',dev_err(k),err(k));

	save([run_name '.mat'],'labels','dev_err','err','C','gamma');
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [f, energ, diagn] = scaled_mh(f, opt, gp, x, y, z)
%SCALED_MH  A scaled Metropolis-Hastings sampling for latent values
%
%  Description
%    [F, ENERG, DIAG] = SCALED_MH(F, OPT, GP, X, Y) takes the
%    current latent values F, options structure OPT, Gaussian
%    process structure GP, inputs X and outputs Y. Samples new
%    latent values and returns also energies ENERG and diagnostics
%    DIAG. The latent values are sampled from their conditional
%    posterior p(f|y,th).
%
%    The latent values are whitened with the prior covariance
%    before the sampling. This reduces the autocorrelation and
%    speeds up the mixing of the sampler. See (Neal, 1998) for
%    details on implementation.
%
%    The options structure should include the following fields:
%      repeat              - the number of MH-steps before
%                            returning a single sample (default 10)
%      sample_latent_scale - the scale for the MH-step (default 0.5)
%
%    OPT = SCALED_MH() Returns default options
%
%    OPT = SCALED_MH(OPT) Returns default options for fields not
%    yet set in OPT
%
%   Reference:
%
%    Neal, R. M. (1998) Regression and classification using
%    Gaussian process priors (with discussion), in J. M. Bernardo,
%    et al (editors) Bayesian Statistics 6, Oxford University
%    Press, pp. 475-501:
%
%  See also
%    GP_MC
%
% Copyright (c) 1999,2011 Aki Vehtari
% Copyright (c) 2006-2010 Jarno Vanhatalo

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

% set default options using hmc2_opt
  if nargin<=1
    if nargin==0
      f=scaled_mh_opt();
    else
      f=scaled_mh_opt(f);
    end
    return
  end

  [n,nout] = size(y);
  if isfield(gp.lik, 'nondiagW')
    switch gp.lik.type
      case {'LGP', 'LGPC'}
        % Do nothing
      case {'Softmax', 'Multinom'}
        % Do nothing
      otherwise
        nout=length(gp.comp_cf);
    end
    if isfield(gp, 'comp_cf')  % own covariance for each ouput component
      multicf = true;
      if length(gp.comp_cf) ~= nout
        error('SCALED_MH: the number of component vectors in gp.comp_cf must be the same as number of outputs.')
      end
    else
      multicf = false;
    end
  end
  f = reshape(f,n,nout);

  maxcut = -log(eps);
  mincut = -log(1/realmin - 1);
  lvs=opt.sample_latent_scale;
  a = max(min(f, maxcut),mincut);

  switch gp.type
    case {'FULL'}

      if ~isfield(gp.lik, 'nondiagW') || ismember(gp.lik.type, {'LGP' 'LGPC'})
        [K,C]=gp_trcov(gp, x);
        if isfield(gp,'meanf')
          [H_m,b_m,B_m]=mean_prep(gp,x,[]);
          C = C + H_m'*B_m*H_m;
        end
        L=chol(C)';
      else
        L = zeros(n,n,nout);
        if multicf
          for i1=1:nout
            [tmp, C] = gp_trcov(gp, x, gp.comp_cf{i1});
            L(:,:,i1)=chol(C, 'lower');
          end
        else
          for i1=1:nout
            [tmp, C] = gp_trcov(gp, x);
            L(:,:,i1)=chol(C, 'lower');
          end
        end
      end
      e = -gp.lik.fh.ll(gp.lik, y, f, z);
      ft = zeros(size(y));

      % Adaptive control algorithm to find such a value for lvs
      % that the rejection rate of Metropolis is optimal.
      slrej = 0;
      for li=1:100
        for i1 =1:nout
          ft(:,i1)=sqrt(1-lvs.^2).*f(:,i1)+lvs.*L(:,:,i1)*randn(n,1);
        end
        ed = -gp.lik.fh.ll(gp.lik, y, ft, z);
        a=e-ed;
        if exp(a) > rand(1)
          f=ft;
          e=ed;
          lvs=min(1,lvs*1.1);
        else
          lvs=max(1e-8,lvs/1.05);
        end
      end
      opt.sample_latent_scale=lvs;
      % Do the actual sampling
      for li=1:(opt.repeat)
        for i1 =1:nout
          ft(:,i1)=sqrt(1-lvs.^2).*f(:,i1)+lvs.*L(:,:,i1)*randn(n,1);
        end
        ed = -gp.lik.fh.ll(gp.lik, y, ft, z);
        a=e-ed;
        if exp(a) > rand(1)
          f=ft;
          e=ed;
        else
          slrej=slrej+1;
        end
      end
      diagn.rej = slrej/opt.repeat;
      diagn.lvs = lvs;
      diagn.opt=opt;
      energ=[];
      f = f(:)';

    case 'FIC'
      u = gp.X_u;
      m = size(u,1);
      % Turn the inducing vector on right direction
      if size(u,2) ~= size(x,2)
        u=u';
      end
      % Calculate some help matrices
      [Kv_ff, Cv_ff] = gp_trvar(gp, x);
      K_fu = gp_cov(gp, x, u);
      K_uu = gp_trcov(gp, u);
      Luu = chol(K_uu)';

      % Evaluate the Lambda (La)
      B=Luu\(K_fu');
      Qv_ff=sum(B.^2)';
      Lav = Cv_ff-Qv_ff;
      sLav = sqrt(Lav);

      n=length(y);
      e = -gp.lik.fh.ll(gp.lik, y, f, z);

      % Adaptive control algorithm to find such a value for lvs
      % so that the rejection rate of Metropolis is optimal.
      slrej = 0;
      for li=1:100
        ft=sqrt(1-lvs.^2).*f + lvs.*(sLav.*randn(n,1) + B'*randn(m,1));
        ed = -gp.lik.fh.ll(gp.lik, y, ft, z);
        a=e-ed;
        if exp(a) > rand(1)
          f=ft;
          e=ed;
          lvs=min(1,lvs*1.1);
        else
          lvs=max(1e-8,lvs/1.05);
        end
      end
      opt.sample_latent_scale=lvs;
      % Do the actual sampling
      for li=1:(opt.repeat)
        ft=sqrt(1-lvs.^2).*f + lvs.*(sLav.*randn(n,1) + B'*randn(m,1));
        ed = -gp.lik.fh.ll(gp.lik, y, ft, z);
        a=e-ed;
        if exp(a) > rand(1)
          f=ft;
          e=ed;
        else
          slrej=slrej+1;
        end
      end
      diagn.rej = slrej/opt.repeat;
      diagn.lvs = lvs;
      diagn.opt=opt;
      energ=[];
      f = f';

    case 'PIC'
      u = gp.X_u;
      m = size(u,1);
      ind = gp.tr_index;
      if size(u,2) ~= size(x,2)
        u=u';
      end

      % Calculate some help matrices
      [Kv_ff, Cv_ff] = gp_trvar(gp, x);  % 1 x f  vector
      K_fu = gp_cov(gp, x, u);         % f x u
      K_uu = gp_trcov(gp, u);    % u x u, noiseles covariance K_uu
      Luu = chol(K_uu)';

      % Evaluate the Lambda (La) for specific model
      % Q_ff = K_fu*inv(K_uu)*K_fu'
      % Here we need only the diag(Q_ff), which is evaluated below
      B=Luu\K_fu';
      iLaKfu = zeros(size(K_fu));  % f x u
      for i=1:length(ind)
        Qbl_ff = B(:,ind{i})'*B(:,ind{i});
        [Kbl_ff, Cbl_ff] = gp_trcov(gp, x(ind{i},:));
        La{i} = Cbl_ff - Qbl_ff;
        CLa{i} = chol(La{i})' ;
      end

      n=length(y);
      e = -gp.lik.fh.ll(gp.lik, y, f, z);

      % Adaptive control algorithm to find such a value for lvs
      % so that the rejection rate of Metropolis is optimal.
      slrej = 0;
      for li=1:100
        sampf = randn(size(f));
        for i=1:length(ind)
          sampf(ind{i},:) = CLa{i}*sampf(ind{i},:);
        end
        ft=sqrt(1-lvs.^2).*f + lvs.*(sampf + B'*randn(m,1));
        at = max(min(ft, maxcut),mincut);
        ed = -gp.lik.fh.ll(gp.lik, y, ft, z);
        a=e-ed;
        if exp(a) > rand(1)
          f=ft;
          e=ed;
          lvs=min(1,lvs*1.1);
        else
          lvs=max(1e-8,lvs/1.05);
        end
      end
      opt.sample_latent_scale=lvs;
      % Do the actual sampling
      for li=1:(opt.repeat)
        sampf = randn(size(f));
        for i=1:length(ind)
          sampf(ind{i},:) = CLa{i}*sampf(ind{i},:);
        end
        ft=sqrt(1-lvs.^2).*f + lvs.*(sampf + B'*randn(m,1));
        at = max(min(ft, maxcut),mincut);
        ed = -gp.lik.fh.ll(gp.lik, y, ft, z);
        a=e-ed;
        if exp(a) > rand(1)
          f=ft;
          e=ed;
        else
          slrej=slrej+1;
        end
      end
      diagn.rej = slrej/opt.repeat;
      diagn.lvs = lvs;
      diagn.opt=opt;
      energ=[];
      f = f';

    case 'CS+FIC'
      u = gp.X_u;
      cf_orig = gp.cf;
      ncf = length(gp.cf);
      n = size(x,1); m = size(u,1);

      cf1 = {};
      cf2 = {};
      j = 1;
      k = 1;
      for i = 1:ncf
        if ~isfield(gp.cf{i},'cs')
          cf1{j} = gp.cf{i};
          j = j + 1;
        else
          cf2{k} = gp.cf{i};
          k = k + 1;
        end
      end
      gp.cf = cf1;

      % First evaluate the needed covariance matrices
      % if they are not in the memory
      % v defines that parameter is a vector
      [Kv_ff, Cv_ff] = gp_trvar(gp, x);  % 1 x f  vector
      K_fu = gp_cov(gp, x, u);         % f x u
      K_uu = gp_trcov(gp, u);          % u x u, noiseles covariance K_uu
      K_uu = (K_uu+K_uu')./2;          % ensure the symmetry of K_uu
      Luu = chol(K_uu)';

      B=Luu\(K_fu');
      Qv_ff=sum(B.^2)';
      Lav = Cv_ff-Qv_ff;   % 1 x f, Vector of diagonal elements
      gp.cf = cf2;
      K_cs = gp_trcov(gp,x);
      La = sparse(1:n,1:n,Lav,n,n) + K_cs;
      gp.cf = cf_orig;

      LD = ldlchol(La);
      sLa = chol(La)';

      n=length(y);
      e = -gp.lik.fh.ll(gp.lik, y, f, z);

      % Adaptive control algorithm to find such a value for lvs
      % so that the rejection rate of Metropolis is optimal.
      slrej = 0;
      for li=1:100
        ft=sqrt(1-lvs.^2).*f + lvs.*(sLa*randn(n,1) + B'*randn(m,1));
        ed = -gp.lik.fh.ll(gp.lik, y, ft, z);
        a=e-ed;
        if exp(a) > rand(1)
          f=ft;
          e=ed;
          lvs=min(1,lvs*1.1);
        else
          lvs=max(1e-8,lvs/1.05);
        end
      end
      opt.sample_latent_scale=lvs;
      % Do the actual sampling
      for li=1:(opt.repeat)
        ft=sqrt(1-lvs.^2).*f + lvs.*(sLa*randn(n,1) + B'*randn(m,1));
        ed = -gp.lik.fh.ll(gp.lik, y, ft, z);
        a=e-ed;
        if exp(a) > rand(1)
          f=ft;
          e=ed;
        else
          slrej=slrej+1;
        end
      end
      diagn.rej = slrej/opt.repeat;
      diagn.lvs = lvs;
      diagn.opt=opt;
      energ=[];
      f = f';

  end
end

function opt = scaled_mh_opt(opt)
%SCALED_MH_OPT  Default options for scaled Metropolis-Hastings sampling
%
%  Description
%    OPT = SCALED_MH_OPT
%      return default options
%    OPT = SCALED_MH_OPT(OPT)
%      fill empty options with default values
%
%  The options and defaults are
%      repeat              - the number of MH-steps before
%                            returning a single sample (default 10)
%      sample_latent_scale - the scale for the MH-step (default 0.5)

  if nargin < 1
    opt=[];
  end

  if ~isfield(opt,'repeat')
    opt.repeat=10;
  end
  if ~isfield(opt,'sample_latent_scale')
    opt.sample_latent_scale=0.5;
  end

end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function Offspring = Operator(Problem,Particle,Pbest,Gbest)
% Particle swarm optimization in SMPSO

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

    %% Parameter setting
    ParticleDec = Particle.decs;
    PbestDec    = Pbest.decs;
    GbestDec    = Gbest.decs;
    [N,D]       = size(ParticleDec);
    ParticleVel = Particle.adds(zeros(N,D));

    %% Particle swarm optimization
    W  = repmat(unifrnd(0.1,0.5,N,1),1,D);
    r1 = repmat(rand(N,1),1,D);
    r2 = repmat(rand(N,1),1,D);
    C1 = repmat(unifrnd(1.5,2.5,N,1),1,D);
    C2 = repmat(unifrnd(1.5,2.5,N,1),1,D);
    OffVel = W.*ParticleVel + C1.*r1.*(PbestDec-ParticleDec) + C2.*r2.*(GbestDec-ParticleDec);
    phi    = max(4,C1+C2);
    OffVel = OffVel.*2./abs(2-phi-sqrt(phi.^2-4*phi));
    delta  = repmat((Problem.upper-Problem.lower)/2,N,1);
    OffVel = max(min(OffVel,delta),-delta);
    OffDec = ParticleDec + OffVel;

    %% Deterministic back
    Lower  = repmat(Problem.lower,N,1);
    Upper  = repmat(Problem.upper,N,1);
    repair = OffDec < Lower | OffDec > Upper;
    OffVel(repair) = 0.001*OffVel(repair);
    OffDec = max(min(OffDec,Upper),Lower);

    %% Polynomial mutation
    disM  = 20;
    Site1 = repmat(rand(N,1)<0.15,1,D);
    Site2 = rand(N,D) < 1/D;
    mu    = rand(N,D);
    temp  = Site1 & Site2 & mu<=0.5;
    OffDec(temp) = OffDec(temp)+(Upper(temp)-Lower(temp)).*((2.*mu(temp)+(1-2.*mu(temp)).*...
                   (1-(OffDec(temp)-Lower(temp))./(Upper(temp)-Lower(temp))).^(disM+1)).^(1/(disM+1))-1);
    temp  = Site1 & Site2 & mu>0.5;
    OffDec(temp) = OffDec(temp)+(Upper(temp)-Lower(temp)).*(1-(2.*(1-mu(temp))+2.*(mu(temp)-0.5).*...
                   (1-(Upper(temp)-OffDec(temp))./(Upper(temp)-Lower(temp))).^(disM+1)).^(1/(disM+1)));
    Offspring = Problem.Evaluation(OffDec,OffVel);
end


text:
function setdefaultbutton(figHandle, btnHandle)
% WARNING: This feature is not supported in MATLAB and the API and
% functionality may change in a future release.

%SETDEFAULTBUTTON Set default button for a figure.
%  SETDEFAULTBUTTON(BTNHANDLE) sets the button passed in to be the default button
%  (the button and callback used when the user hits "enter" or "return"
%  when in a dialog box.
%
%  This function is used by inputdlg.m, msgbox.m, questdlg.m and
%  uigetpref.m.
%
%  Example:
%
%  f = figure;
%  b1 = uicontrol('style', 'pushbutton', 'string', 'first', ...
%       'position', [100 100 50 20]);
%  b2 = uicontrol('style', 'pushbutton', 'string', 'second', ...
%       'position', [200 100 50 20]);
%  b3 = uicontrol('style', 'pushbutton', 'string', 'third', ...
%       'position', [300 100 50 20]);
%  setdefaultbutton(b2);
%

%  Copyright 2005-2007 The MathWorks, Inc.

%--------------------------------------- NOTE ------------------------------------------
% This file was copied into matlab/toolbox/local/private.
% These two files should be kept in sync - when editing please make sure
% that *both* files are modified.

% Nargin Check
if nargin<1, error('MATLAB:setdefaultbutton:InvalidNumberOfArguments','Too few arguments for setdefaultbutton'); end
if nargin>2, error('MATLAB:setdefaultbutton:InvalidNumberOfArguments','Too many arguments for setdefaultbutton'); end

try
    ujava = feature('javafigures');
catch
    ujava = true;
end
if ujava && (usejava('awt') == 1)
    % We are running with Java Figures
    useJavaDefaultButton(figHandle, btnHandle)
else
    % We are running with Native Figures
    useHGDefaultButton(figHandle, btnHandle);
end

    function useJavaDefaultButton(figH, btnH)
        % Get a UDD handle for the figure.
        fh = handle(figH);
        % Call the setDefaultButton method on the figure handle
        fh.setDefaultButton(btnH);
    end

    function useHGDefaultButton(figHandle, btnHandle)
        % First get the position of the button.
        btnPos = getpixelposition(btnHandle);

        % Next calculate offsets.
        leftOffset   = btnPos(1) - 1;
        bottomOffset = btnPos(2) - 2;
        widthOffset  = btnPos(3) + 3;
        heightOffset = btnPos(4) + 3;

        % Create the default button look with a uipanel.
        % Use black border color even on Mac or Windows-XP (XP scheme) since
        % this is in natve figures which uses the Win2K style buttons on Windows
        % and Motif buttons on the Mac.
        h1 = uipanel(get(btnHandle, 'Parent'), 'HighlightColor', 'black', ...
            'BorderType', 'etchedout', 'units', 'pixels', ...
            'Position', [leftOffset bottomOffset widthOffset heightOffset]);

        % Make sure it is stacked on the bottom.
        uistack(h1, 'bottom');
    end
end

author: spm


text:
function t_obj = parcel_stats2statistic_image(atlas_obj, tscores, pvalues, dfe, sig)
% Utility that transforms parcel-wise t-statistics and p-values into a statistic_image object in voxelwise image space
% Expands parcel values so that they are replicated for each voxel in the parcel.
%
% :Usage:
% ::
%
% t_obj = parcel_stats2statistic_image(atlas_obj, tscores, pvalues, dfe, sig)
%
% For objects: Type methods(object_name) for a list of special commands
%              Type help object_name.method_name for help on specific
%              methods.
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2021 Tor Wager
%
%     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/>.
% ..
%
% :Inputs:
%
%   **atlas_obj:**
%        An atlas-class object with k distinct regions or parcels
%
%   **tscores:**
%        n x k matrix of k t-scores for each of n parcels, where k is,
%        e.g., the number of predictors tested in a multiple regression
%
%   **pvalues:**
%        n x k matrix of k p-values for each of n parcels
%
%   **dfe:**
%        n x 1 matrix of error degrees of freedom for each of n parcels
%
%   **sig:**
%        n x k matrix of k significant tests for each of n parcels
%        This may be, e.g., corrected for multiple comparisons.
%
% :Optional Inputs:
%
% :Outputs:
%
%   **t_obj:**
%        A statistic_image object containing t-maps. .dat contains
%        t-scores, .p contains p-values. Use statistic_image methods like
%        montage(t_obj), orthviews(t_obj), threshold(t_obj) to interact
%        with this object.
%
%
% :Examples:
% ::
%
% :References:
%   None
%
% :See also:
%   parcel_data2fmri_data
%

% Create placeholder statistic image

atlas_obj = replace_empty(atlas_obj);
k = size(tscores, 2);
placeholder_vec = ones(atlas_obj.volInfo.n_inmask, k);
all_parcel_idx = double(atlas_obj.dat);
u = unique(all_parcel_idx); u(u == 0) = [];

% initialize variables with correct size

t_obj = statistic_image('dat', single(0 .* placeholder_vec), ...
    'p', placeholder_vec, ...
    'sig', logical(placeholder_vec), ...
    'type', 'T', ...
    'dfe', placeholder_vec, ...
    'volInfo', atlas_obj.volInfo);

% number of parcels

n = size(tscores, 1); % num parcels

if n ~= length(u)
        error('Parcel indices in atlas_obj and extracted parcel-wise t-values do not match. Check code.')
end

t_obj.dat = zeros(size(placeholder_vec, 1), k); % voxels x conditions

    for j = 1:length(u)
        % For each parcel, fill in statistic_image object
        % ----------------------------------------------------
        parcelidx = u(j);

        wh_vox = all_parcel_idx == parcelidx;

        % map parcels to voxels
        t_obj.dat(wh_vox, :) = repmat(tscores(j, :), sum(wh_vox), 1);
        t_obj.p(wh_vox, :) = repmat(pvalues(j, :), sum(wh_vox), 1);
        t_obj.dfe(wh_vox, :) = repmat(dfe(j, 1), sum(wh_vox), k);       % replicate dfe
        t_obj.sig(wh_vox, :) = repmat(sig(j, :), sum(wh_vox), 1);

    end

t_obj = enforce_variable_types(t_obj);  % space-saving: 5/24/17

% input_struct.t_statistic_obj = t_obj;

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [Ht,weights] = riskmetrics2006(data,tau0,tau1,kmax,rho)
% Computes the Riskmetrics 2006 covariance, which is a weighted average of
% EWMA covariances
%
% USAGE:
%  [HT,WEIGHTS] = riskmetrics2006(DATA,TAU0,TAU1,KMAX,RHO)
%
% INPUTS:
%   DATA     - A T by K matrix of zero mean residuals -OR-
%                K by K by T array of covariance estimators (e.g. realized covariance)
%   TAU0     - [OPTIONAL] Half-life of slowest EWMA
%   TAU1     - [OPTIONAL] Half-life of fastest EWMA
%   KMAX     - [OPTIONAL] Number of EWMA components to use
%   RHO      - [OPTIONAL] Decay factor to use in half-lives.  The
%                half-lives used are TAU1, TAU1*RHO, TAU1*RHO^2, ...
%
% OUTPUTS:
%   HT       - A [K K T] dimension matrix of conditional covariances
%   WEIGHTS  - A T by T matrix which contains the final weights used
%                computing all conditional covariances.  WEIGHT(i,j)
%                contains the weight of the outer-product in time period j
%                in the covariance forecast at period i+1
%
% COMMENTS:
%   The conditional variance, H(t), of a RM2006 covariance model
%
%      H(t) = w(i)*Htilde(t,i)
%
%   where Htilde(t,i) is an EWMA covariance, i=1,2,..., and
%
%   w(i) = 1-log(TAU1^((i-1)*RHO))/log(TAU0)
%
%   where the weights have been normalized so that they sum to 1.  All
%   EWMAs are initialized using a backward EWMA with the same decay.
%
% EXAMPLES:
%   RiskMetrics 2006 methodology using the suggested reference values
%     tau0   = 1560;
%     tau1   = 4;
%     taumax = 512;
%     rho    = sqrt(2);
%     kmax   = round(log(taumax/tau1)/log(rho)); % 14
%     Ht     = riskmetrics2006(data,tau0,tau1,kmax,rho)
%
%   RiskMetrics 1994 methodology using .94 as a special case of RM2006
%     tau0 = 1560;         % Does not matter
%     rho  = 1;            % Does not matter
%     tau1 = -1/log(.94);
%     kmax = 1;
%     Ht = riskmetrics2006(data,tau0,tau1,kmax,rho)

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 03/10/2011

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch nargin
    case 1
        tau0=[];
        tau1=[];
        kmax=[];
        rho =[];
    case 2
        tau1=[];
        kmax=[];
        rho =[];
    case 3
        kmax=[];
        rho =[];
    case 4
        rho =[];
    case 5
        % Nothing
    otherwise
        error('1 to 5 inputs required.')
end

if isempty(tau0)
    tau0 = 1560;
end
if isempty(tau1)
    tau1 = 4;
end
if isempty(kmax)
    kmax = 14;
end
if isempty(rho)
    rho = sqrt(2);
end

if tau1*rho^(kmax-1)>tau0
    error('The inputs must satisfy: TAU1*RHO^(KMAX-1)<TAU0')
end
if tau1<0
    error('TAU1 must be positive')
end
if tau0<0
    error('TAU0 must be positive')
end
if kmax<1 || floor(kmax)~=kmax
    error('KMAX must be an integer (weakly) larger than 1.')
end
if rho<0
    error('RHO must be positive')
end

if ndims(data)==2
    [T,K] = size(data);
    temp = zeros(K,K,T);
    for t=1:T
        temp(:,:,t) = data(t,:)'*data(t,:);
    end
    data = temp;
else
    K = size(data,1);
    T = size(data,3);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tauks = tau1*rho.^((1:kmax)-1);
w= 1-log(tauks)/log(tau0);
w=w/sum(w);
Ht = zeros(K,K,T);
Httilde = zeros(K,K,T);

for k=1:kmax
    tauk = tauks(k);
    mu = exp(-1/tauk);
    % back casting
    endPoint = max(min(floor(log(.01)/log(mu)),T),k);
    weights = (1-mu).*mu.^(0:endPoint-1);
    weights = weights/sum(weights);
    backCast = zeros(K);
    for i=1:endPoint
        backCast = backCast + weights(i)*data(:,:,i);
    end
    Httilde(:,:,1) = backCast;

    for t=2:T
        Httilde(:,:,t) = mu*Httilde(:,:,t-1) + (1-mu)*data(:,:,t-1);
    end
    Ht = Ht + w(k) * Httilde;
end

if nargout>1
    weights = zeros(T,T);
    for k=1:kmax
        tauk = tauks(k);
        mu = exp(-1/tauk);
        weightMatrix = [mu.^(0:T-1)' zeros(T,T-1)];
        weightMatrix(T,T) = (1-mu);
        for j=1:(T-2);
            weightMatrix(T,T-j) = mu * weightMatrix(T,T-j+1);
        end
        for t=2:T-1
            weightMatrix(t,2:t) = weightMatrix(T,T+(-t+2:0));
        end
        weights = weights + w(k) * weightMatrix;
    end
end

author: bashtage

repo: mfe-toolbox


text:
classdef AbstractFixedStepIntegrator < AbstractIntegrator
    %AbstractFixedStepIntegrator Summary of this class goes here
    %   Detailed explanation goes here

    properties(Constant, Access=private)
        %RK5
%         At = [ 1/5,          0,           0,            0,         0
%                3/40,         9/40,        0,            0,         0
%                44/45        -56/15,       32/9,         0,         0
%                19372/6561,  -25360/2187,  64448/6561,  -212/729,   0
%                9017/3168,   -355/33,      46732/5247,   49/176,   -5103/18656];
%
%         Bt = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84];
%
%         Ct = [1/5; 3/10; 4/5; 8/9; 1];

        %RK4
%         At = [1/2,  0,  0;
%                0, 1/2,  0;
%                0,   0,  1;]
%
%        Bt = [1/6, 1/3, 1/3, 1/6];
%
%        Ct = [1/2, 1/2, 1];

        %RK2
%         At = [1, 1];
%         Bt = [1/2, 1/2];
%
%         Ct = [1];

        %Ralston's 4th Order Method
%         At = [1/2,   0,   0;
%                 0, 1/2,   0;
%                 0,   0,   1];
%
%        Bt = [1/6, 1/3, 1/3, 1/6];
%
%        Ct = [1/2, 1/2, 1];
    end

    methods(Static)
        function [t,y, te,ye,ie] = integrate(odefun, tspan, y0, options)
            [h,y0,~] = AbstractFixedStepIntegrator.checkInputs(odefun, tspan, y0);

            if(all(h>0))
                propDirPos = true;
            elseif(all(h<0))
                propDirPos = false;
            else
                error('Prop direction might be consistently forward or backwards.');
            end

            if(not(isempty(options.OutputFcn)) && ...
               isa(options.OutputFcn,'function_handle'))
                hasOutput = true;
                outputFcn = options.OutputFcn;
            else
                hasOutput = false;
                outputFcn = [];
            end

            if(not(isempty(options.Events)) && ...
               isa(options.Events,'function_handle'))
                hasEvents = true;
                eventFcn = options.Events;
            else
                hasEvents = false;
                eventFcn = [];
            end

            neq = length(y0);
            N = length(tspan);
            t = NaN(N,1);
            Y = NaN(neq,N);

            Y(:,1) = y0;
            t(1) = tspan(1);

            if(hasOutput)
                outputFcn([tspan(1) tspan(end)],y0,'init');
            end

            ie = [];
            te = [];
            ye = [];
            if(hasEvents)
                [prevEvtValues,~,~] = eventFcn(t(1),Y(:,1)); %#ok<RHSFN>
            else
                prevEvtValues = [];
            end

            evtTerminate = false;
            outputTerminate = false;

            i = 2;
            while(keepLooping(propDirPos, i, tspan))
                ti = tspan(i-1);
                if(breakLoop(propDirPos, ti, tspan(end)))
                    break;
                end

                hi = h(i-1);
                yi = Y(:,i-1);

                t(i) = tspan(i);
                Y(:,i) = AbstractFixedStepIntegrator.stepOnce(odefun, ti, yi, hi, neq);

                if(hasEvents)
                    [evtValues,evtIsTerminal,evtDirection] = eventFcn(t(i),Y(:,i)); %#ok<RHSFN>

                    signChangesBool = (prevEvtValues .* evtValues < 0);
                    evtValZeroBool = evtValues == 0;

                    rootDetect = signChangesBool | evtValZeroBool;
                    if(any(rootDetect))
                        rootFcnFzero = @(subHi) AbstractFixedStepIntegrator.rootFindingFunc(subHi, odefun, ti, yi, neq, eventFcn, true);

                        if(all(rootDetect == evtValZeroBool)) %we nailed it on the money
                            tSubI = t(i);
                            exitflag = 1;

                        else %use root finding to actually find the event root
                            m = (evtValues - prevEvtValues) / (t(i) - t(i-1));
                            b = (evtValues - t(i) .* m);
                            t0 = -b./m;

                            t0(rootDetect == 0) = NaN;
                            [minT0, ~] = min(t0);

%                             [tSubI, ~, exitflag] = fzero(rootFcnFzero, minT0);

                            tSubI = SecantMethod(rootFcnFzero, t(i-1), minT0, 1e-6);

                            if(isnan(tSubI))
                                exitflag = -1;
                            else
                                exitflag = 1;
                            end

                            if(exitflag ~= 1 || tSubI < t(i-1) || tSubI > t(i))
                                rootFcnFminBnd = @(subHi) AbstractFixedStepIntegrator.rootFindingFunc(subHi, odefun, ti, yi, neq, eventFcn, false);
                                [tSubI,~,exitflag,~] = fminbnd(rootFcnFminBnd, t(i-1), t(i), optimset('TolX',1E-5));
                            end
                        end

                        if(exitflag == 1)
                            diffForDeriv = sign(evtValues - prevEvtValues);

                            [~, evtInd] = rootFcnFzero(tSubI);
                            if(evtDirection(evtInd) == 0 || ...
                               diffForDeriv(evtInd) == evtDirection(evtInd))

                                tspan1 = tspan(1:i-1);
                                tspan2 = tspan(i:end);
                                tspan = [tspan1, tSubI, tspan2];

                                h1 = h(1:i-1);
                                h2 = h(i:end);
                                hi = tSubI - ti;
                                h = [h1, hi, h2];

                                t(i) = tspan(i);
                                Y(:,i) = AbstractFixedStepIntegrator.stepOnce(odefun, ti, yi, hi, neq);

                                ie(end+1,1) = evtInd;
                                te(end+1,1) = tSubI;
                                ye(end+1,:) = Y(:,i)';

                                if(evtIsTerminal(evtInd) == 1)
                                    evtTerminate = true;
                                else
                                    evtTerminate = false;
                                end
                            else
                                evtTerminate = false;
                            end
                        end
                    end

                    prevEvtValues = evtValues;
                end

                if(hasOutput)
                    status = outputFcn(t(i),Y(:,i),[]);
                    if(status == 1)
                        outputTerminate = true;
                    else
                        outputTerminate = false;
                    end
                end

                if(evtTerminate || outputTerminate)
                    break;
                end

                i=i+1;
            end

            y = Y.';

            bool = ~isnan(t);
            t = t(bool);
            y = y(bool,:);

            if(hasOutput)
                outputFcn([],[],'done');
            end
        end

    end

    methods(Static, Access=protected)
        [A,B,C] = getButcherTableauData();

        function [h,y0,f0] = checkInputs(odefun, tspan, y0)
            if ~isnumeric(tspan)
                error('TSPAN should be a vector of integration steps.');
            end

            if ~isnumeric(y0)
                error('Y0 should be a vector of initial conditions.');
            end

            h = diff(tspan);
            if any(sign(h(1))*h <= 0)
                error('Entries of TSPAN are not in order.')
            end

            y0 = y0(:);   % Make a column vector.
            try
                f0 = odefun(tspan(1),y0);
            catch ME
                msg = ['Unable to evaluate the ODEFUN at t0,y0: ',ME.message];
                error(msg);
            end

            if ~isequal(size(y0),size(f0))
                error('Inconsistent sizes of Y0 and f(t0,y0).');
            end
        end

        function [yNp1, fEvals] = stepOnce(odefun, ti, yi, hi, neq)
            [A,B,C] = ODE5Integrator.getButcherTableauData(); %needs to be generalized

            nstages = length(B);
            F = NaN(neq,nstages);

            % General explicit Runge-Kutta framework
            F(:,1) = odefun(ti,yi);
            for stage = 2:nstages
                tstage = ti + C(stage-1)*hi;
                ystage = yi + F(:,1:stage-1)*(hi*A(1:stage-1,stage-1));
                F(:,stage) = odefun(tstage,ystage);
            end
            yNp1 = yi + F*(hi*B);
            fEvals = 1 + (nstages - 1);
        end
    end

    methods(Static, Access=private)
        function [rootOutput, I] = rootFindingFunc(subT, odefun, ti, yi, neq, eventFcn, forRootFinder)
            subHi = subT - ti;
            yNp1 = AbstractFixedStepIntegrator.stepOnce(odefun, ti, yi, subHi, neq);

            [evtValues,~,~] = eventFcn(subT,yNp1);

            rootOutput = evtValues;
            [~,I] = min(abs(rootOutput));

            if(forRootFinder)
                rootOutput = rootOutput(I); %for fzero and the like
            else
                rootOutput = abs(rootOutput(I)); %for an optimizer like fminbnd
            end
        end
    end
end

function tf = keepLooping(propDirPos, i, tspan)
    ti = tspan(i-1);
    if(propDirPos)
        tf = ti < tspan(end);
    else
        tf = ti > tspan(end);
    end
end

function tf = breakLoop(propDirPos, ti, tspanEnd)
    if(propDirPos)
        tf = ti >= tspanEnd;
    else
        tf = ti <= tspanEnd;
    end
end


text:
%kvhmed 'Median Filter the Image Using Histogram to Find Median (K1)'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vhmed.pane file
%
% Parameters:
% InputFile: i 'Input Image ', required: 'input image filename'
% Integer: w 'Window Width ', default: 3: 'window width argument'
% Integer: h 'Window Height', default: 3: 'window height argument'
% OutputFile: o 'Output Image', required: 'output image filename'
%
% Example: o = kvhmed(i, {'i','';'w',3;'h',3;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vhmed - Median Filter the Image Using Histogram to Find Median  (K1)
%
%  DESCRIPTION
% .I vhmed
% computes a two dimensional median filter of width w and height h over
% the given image.  The median is computed by moving a w by h
% pixel mask over the image and updating the histogram of the mask.
%
% Images must be of data storage type byte.
%
% In rectangular odd windows (masks) the median value is replaced in
% the center pixel location of the window. The default window
% size is 3x3.
%
% The width or height of a window cannot be an even number.
%
%
%
%  EXAMPLES
% vhmed -i input.image -o output.image -w 3 -h 3
%
% This command performs median filtering of width 3 and height 3 on
% input.image and stores the result in output.image.
%
%  "SEE ALSO"
% vqmed(1)
%
%  RESTRICTIONS
% .I vhmed
% does not support even window dimensions and only works on
% data storage type byte.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

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

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [data,units] = compute_min_dwing_angle_out(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  data{i} = min(-diff(trx(fly).wing_anglel),diff(trx(fly).wing_angler)) ./ trx(fly).dt;

end
units = parseunits('rad/s');

author: kristinbranson

repo: JAABA


text:
function [F, subF, Data] = ObjFun(Track, UAV)
%OBJFUN 目标函数、适应度函数（一个agent的）

% 多目标优化权重设置（必需为行向量）
weight = [ 0.05, 0.05, 0.1, 0.7, 0.7 ]; % 默认权重

% 表达式系数调整（将各项指标无量纲化）
p1 = 1;   % 燃料项（已除以最大航程）
p21 = 1; % 高度项
p22 = 1; % 低度项
p31 = 1.2; % 雷达威胁
p32 = 1.1; % 其余威胁
p4 = 1.2;   % 时间同步项
p5 = 1;   % 碰撞项

% 进行航迹检测
report = TrackDetect(Track, UAV);  % Track 为 struck结构

% 燃油（论文3）
ZZ = sum(UAV.limt.L);
MaxL_mt = ZZ(2);
f_o = p1 * report.L_mt / MaxL_mt;

% 高度（论文3）
dim = UAV.PointDim;
if dim < 3
    f_h = 0;
else
    f_h = 0;
    for i = 1 : UAV.num
        Hmax = UAV.limt.h(i, 2);
        Hmin = UAV.limt.h(i, 1);
        for k = 1 : UAV.PointNum(i)
            z = Track.P{i}(3, k);
            if z>Hmax
                fk = p21 * (z - Hmax);
            elseif z >= Hmin
                fk = 0;
            else
                fk = p22 * (Hmin - z);
            end
            f_h = f_h + fk;
        end
    end
end

% 威胁（论文2）
% 注意：只适应球或圆区域，不适合圆柱区域
O_r = UAV.Menace.radar(: ,1:end-1);             %雷达
O_o = UAV.Menace.other(: ,1:end-1);            %导弹，火炮，气象等
f_t = 0;  % 威胁代价
for i = 1 : UAV.num
    for k = 1 : UAV.PointNum(i)
        P = Track.P{i}(:, k)' ;  % 转置成 1*dim
        for m = 1 : size(O_r, 1)
            fk = p31 / (norm(P - O_r(m, :)))^4;
            f_t = f_t + fk;
        end
        for m = 1 : size(O_o, 1)
            fk = p32 / norm(P - O_o(m, :));
            f_t = f_t + fk;
        end
    end
end

% 同步（论文1）
f_m = 0;
for i = 1 : UAV.num
    Li = report.L(i);
    tmax = Li / UAV.limt.v(i,1);
    tmin = Li / UAV.limt.v(i,2);
    ti = report.time(i);
    tc = UAV.tc;
    if tc <= tmax && tc >= tmin
        fk = 0;
    else
        fk = p4 * abs(ti - tc);
    end
    f_m = f_m + fk;
end

% 碰撞（论文1）
f_c = p5 * report.col_times;

% 目标函数分量（必需是列向量，否则无法聚类）
subF = [ f_o; f_h; f_t; f_m; f_c ]; % 5*1

% 加权目标函数
F = weight * subF ;

% 输出信息
Data.ProbPoint = report.ProbPoint;      % 所有有问题的点
Data.AngleProb = report.AngleProb;    % 不满足角度约束的点
Data.TrajProb = report.TrajProb;           % 不满足最小航迹间隔的点
Data.Threat = report.Threat;                  % 受威胁的点

Data.L = report.L;                                      % 每个无人机的航程
Data.t = report.time;                                 % 每个无人机的时间
Data.c = report.col_times;                         % 所有无人机总碰撞次数

end

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


text:
function Pout = filename_get_new_root_dir(P,newroot,nlevels)
% :Usage:
% ::
%
%     Pout = filename_get_new_root_dir(P,newroot,nlevels)
%
% :Inputs:
%
%   **P:**
%        is files (string matrix)
%
%   **newroot:**
%        Append a new root directory to filenames
%        if newroot is empty, prompt
%
%   **nlevels:**
%        Keeps nlevels dirs deep from old directory:
%          - keeps only filename
%          - keeps filename plus last subdir
%          - keeps filename plus last 2 subdirs, etc.
%
% :Examples:
% ::
%
%    Pout = filename_get_new_root_dir(EXPT.SNPM.P{3},pwd,2)
%
%    % Change root dir to pwd for all images in EXPT
%    for i=1:length(EXPT.SNPM.P), EXPT.SNPM.P{i} = filename_get_new_root_dir(EXPT.SNPM.P{i},pwd,2); end
%
% Note: see regexprep.m for a simpler way!!! This function could be
% improved by using it.
%
% ..
%    tor wager
% ..

if isempty(newroot), newroot = spm_get(-1,'*','Get new root directory'); end

Pout = get_Pout(P(1,:),newroot,nlevels);

for i = 2:size(P,1)

    Pout = str2mat(Pout,get_Pout(P(i,:),newroot,nlevels));

end

for i = 1:size(Pout, 1)
    isfile(i) = exist(deblank(Pout(i,:)), 'file');
end
fprintf(1,'Found %3.0f valid files.\n', sum(isfile > 0));

return

function Pout = get_Pout(p,str,nlevels);

p = deblank(p);
[d,f,e] = fileparts(p);
filename = [f e];

for i = 2:nlevels
    [d,dd] = fileparts(d);
    filename = fullfile(dd,filename);
end

Pout = fullfile(str,filename);

return


text:
function [new_nbrs, new_ops, new_nodes, new_topos] = mk_nbrs_of_dag_topo(G0)
% MK_NBRS_OF_DAG_TOPO Make all DAGs that differ from G0 by a single edge deletion, addition or reversal
% [new_nbrs, new_ops, new_nodes, new_topos] = mk_nbrs_of_dag_topo(G0)
%
% new_nbrs{i} is the i'th neighbor of G0.
% new_ops{i} = 'add', 'del', or 'rev' is the operation used to create the i'th neighbor.
% new_nodes(i,1:2) are the head and tail of the operated-on arc.
% new_topos are topological orders of the neighbours.
%
% We implement the fast acyclicity check described by P. Giudici and R. Castelo,
% "Improving MCMC model search for data mining", submitted to J. Machine Learning, 2001.
%
% Written by Qian Diao <qian.diao@intel.com> on 19 Nov 01
% Reference are ..\BNT\graph\mk_nbrs_of_dag.m, ..\BNT\learning\learn_struct_mcmc.m and
% ..\BNT\graph\topological_sort.m
% Copyright Intel 2001
%
new_nbrs = {};
new_ops = {};
new_nodes = [];
new_topos = {};
cs = {};

n = length(G0);
indeg = zeros(1,n);
zero_indeg = []; % a stack of nodes with no parents
for i=1:n
  indeg(i) = length(parents(G0,i));
  cs{i} = children(G0, i);
  if indeg(i)==0
    zero_indeg = [i zero_indeg];
  end
end

dag = G0;
[nbrs, ops, nodes] = mk_nbrs_of_digraph(dag);
A = init_ancestor_matrix(dag);
%assert(acyclic(new_dag));

d1 = 1;
for d = 1:length(ops)
  i = nodes(d, 1); j = nodes(d, 2);
  legal = 0;
  switch ops{d}
   case 'add',
    if A(i,j)==0
      legal = 1;
    end
   case 'del',
    legal = 1;
   case 'rev',
    ps = mysetdiff(parents(dag, j), i);
    % if any(A(ps,i)) then there is a path i -> parent of j -> j
    % so reversing i->j would create a cycle
    legal = ~any(A(ps, i));
  end

  if legal
    tmp = nbrs(:,:,d);
    new_nbrs{d1} = tmp;
    new_ops{d1} =  ops{d};
    new_nodes(d1,1:2) = nodes(d,1:2);

    % obtain the topological orders of neighbour dags
    zero_indeg_nbr = [];
    indeg_nbr = [];
    cs_nbr = [];

    switch ops{d}
     case 'add' % i is a new parent of j
      zero_indeg_nbr = zero_indeg;
      indeg_nbr = indeg;
      if ~isempty(find(zero_indeg == j)) % j is not a root anymore
        zero_indeg_nbr = mysetdiff(zero_indeg, j);
      end
      indeg_nbr(j) = indeg(j)+1;

      t_nbr=1;
      order_nbr = zeros(1,n);
      while ~isempty(zero_indeg_nbr)
        v_nbr = zero_indeg_nbr(1); % pop v
        zero_indeg_nbr = zero_indeg_nbr(2:end);
        order_nbr(t_nbr) = v_nbr;
        t_nbr = t_nbr + 1;
        if v_nbr == i % j is a new child of i
          cs_nbr = sort([j cs{i}]);
        else
          cs_nbr = cs{v_nbr};
        end
        for k = 1:length(cs_nbr)
          c_nbr = cs_nbr(k);
          indeg_nbr(c_nbr) = indeg_nbr(c_nbr) - 1;
          if indeg_nbr(c_nbr) == 0
            zero_indeg_nbr = [c_nbr zero_indeg_nbr]; % push c
          end
        end
      end

     case 'del' % i is not a parent of j anymore
      zero_indeg_nbr = zero_indeg;
      indeg_nbr = indeg;
      if length(parents(tmp, j))==0
        zero_indeg_nbr = -sort(-[zero_indeg, j]); % descending order
      end
      indeg_nbr(j) = indeg(j) - 1;

      t_nbr=1;
      order_nbr = zeros(1,n);
      while ~isempty(zero_indeg_nbr)
        v_nbr = zero_indeg_nbr(1); % pop v
        zero_indeg_nbr = zero_indeg_nbr(2:end);
        order_nbr(t_nbr) = v_nbr;
        t_nbr = t_nbr + 1;
        if v_nbr == i % j is not a child of i anymore
          cs_nbr = mysetdiff(cs{i}, j);
        else
          cs_nbr = cs{v_nbr};
        end
        for k = 1:length(cs_nbr)
          c_nbr = cs_nbr(k);
          indeg_nbr(c_nbr) = indeg_nbr(c_nbr) - 1;
          if indeg_nbr(c_nbr) == 0
            zero_indeg_nbr = [c_nbr zero_indeg_nbr]; % push c
          end
        end
      end

     case 'rev' %i is a new child of j and j is a new parent of i
      zero_indeg_nbr = zero_indeg;
      indeg_nbr = indeg;
      if ~isempty(find(zero_indeg == i))
        zero_indeg_nbr = mysetdiff(zero_indeg_nbr, i);
      end
      if length(parents(tmp, j))==0
        zero_indeg_nbr = -sort(-[zero_indeg_nbr, j]); % decending order
      end
      indeg_nbr(i) = indeg(i)+1;
      indeg_nbr(j) = indeg(j)-1;

      t_nbr=1;
      order_nbr = zeros(1,n);
      while ~isempty(zero_indeg_nbr)
        v_nbr = zero_indeg_nbr(1); % pop v
        zero_indeg_nbr = zero_indeg_nbr(2:end);
        order_nbr(t_nbr) = v_nbr;
        t_nbr = t_nbr + 1;
        cs_nbr = cs{v_nbr};
        if v_nbr == i % j is not a child of i anymore
          cs_nbr = mysetdiff(cs{i}, j);
        end
        if v_nbr == j % i is a new child of j
          cs_nbr = sort([i cs{j}]);
        end
        for k = 1:length(cs_nbr)
          c_nbr = cs_nbr(k);
          indeg_nbr(c_nbr) = indeg_nbr(c_nbr) - 1;
          if indeg_nbr(c_nbr) == 0
            zero_indeg_nbr = [c_nbr zero_indeg_nbr]; % push c
          end
        end
      end
     end

    new_topos{d1} = order_nbr;
    d1 = d1+1;
  end
end

clear nbrs ops nodes;

%%%%%%%%%
function A = update_row(A, j, dag)
% We compute row j of A
A(j, :) = 0;
ps = parents(dag, j);
if ~isempty(ps)
  A(j, ps) = 1;
end
for k=ps(:)'
  anck = find(A(k,:));
  if ~isempty(anck)
    A(j, anck) = 1;
  end
end

%%%%%%%%
function A = init_ancestor_matrix(dag)
order = topological_sort(dag);
A = zeros(length(dag));
for j=order(:)'
  A = update_row(A, j, dag);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% pop_erpparams() - Set plotting and statistics parameters for cluster ERP
%                   plotting
% Usage:
%   >> STUDY = pop_erpparams(STUDY, 'key', 'val');
%
% Inputs:
%   STUDY        - EEGLAB STUDY set
%
% Input:
%   'topotime'   - [real] Plot ERP scalp maps at one specific latency (ms).
%                   A latency range [min max] may also be defined (the
%                   ERP is then averaged over the interval) {default: []}
%   'filter'     - [real] low pass filter the ERP curves at a given
%                  frequency threshold. Default is no filtering.
%   'timerange'  - [min max] ERP plotting latency range in ms.
%                  {default: the whole epoch}
%   'ylim'       - [min max] ERP limits in microvolts {default: from data}
%   'plotgroups' - ['together'|'apart'] 'together' -> plot subject groups
%                  on the same axis in different colors, else ('apart')
%                  on different axes. {default: 'apart'}
%   'plotconditions' - ['together'|'apart'] 'together' -> plot conditions
%                  on the same axis in different colors, else ('apart')
%                  on different axes. Note: Keywords 'plotgroups' and
%                  'plotconditions' may not both be set to 'together'.
%                  {default: 'together'}
%   'averagechan' - ['on'|'off'] average data channels when several are
%                  selected.
%
% See also: std_erpplot()
%
% Authors: Arnaud Delorme, CERCO, CNRS, 2006-

% Copyright (C) Arnaud Delorme, 2006
%
% 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

function [ STUDY, com ] = pop_erpparams(STUDY, varargin);

STUDY = default_params(STUDY);
TMPSTUDY = STUDY;
com = '';
if isempty(varargin)

    enablecond         = fastif(length(STUDY.design(STUDY.currentdesign).variable(1).value)>1, 'on', 'off');
    enablegroup        = fastif(length(STUDY.design(STUDY.currentdesign).variable(2).value)>1, 'on', 'off');
    detachplots        = fastif(strcmpi(STUDY.etc.erpparams.detachplots,'on'), 1, 0);
    plotconditions     = fastif(strcmpi(STUDY.etc.erpparams.plotconditions, 'together'), 1, 0);
    plotgroups         = fastif(strcmpi(STUDY.etc.erpparams.plotgroups,'together'), 1, 0);
    radio_averagechan  = fastif(strcmpi(STUDY.etc.erpparams.averagechan,'on'), 1, 0);
    radio_scalptopo    = fastif(isempty(STUDY.etc.erpparams.topotime), 0, 1);
    if radio_scalptopo, radio_averagechan = 0; end;
    if radio_scalptopo+radio_averagechan == 0, radio_scalparray = 1; else radio_scalparray = 0; end;

    cb_radio = [ 'set(findobj(gcbf, ''userdata'', ''radio''), ''value'', 0);' ...
                 'set(gcbo, ''value'', 1);' ...
                 'set(findobj(gcbf, ''tag'', ''topotime''), ''string'', '''');' ];
    cb_edit  = [ 'set(findobj(gcbf, ''userdata'', ''radio''), ''value'', 0);' ...
                 'set(findobj(gcbf, ''tag'', ''scalptopotext''), ''value'', 1);' ];
    uilist = { ...
           {'style' 'text'       'string' 'ERP plotting options' 'fontweight' 'bold' 'fontsize', 12} ...
        {} {'style' 'text'       'string' 'Time limits (ms) [low high]' } ...
           {'style' 'edit'       'string' num2str(STUDY.etc.erpparams.timerange) 'tag' 'timerange' } ...
        {} {'style' 'text'       'string' 'Plot limits [low high]'} ...
           {'style' 'edit'       'string' num2str(STUDY.etc.erpparams.ylim) 'tag' 'ylim' } ...
        {} {'style' 'text'       'string' 'Lowpass plotted data [Hz]' } ...
           {'style' 'edit'       'string' num2str(STUDY.etc.erpparams.filter) 'tag' 'filter' } ...
        {} ...
           {'style' 'text'       'string' 'ERP plotting format' 'fontweight' 'bold' 'fontsize', 12} ...
        {} {'style' 'checkbox'   'string' 'Plot first variable on the same panel' 'value' plotconditions 'enable' enablecond  'tag' 'plotconditions' } ...
        {} {'style' 'checkbox'   'string' 'Plot second variable on the same panel' 'value' plotgroups 'enable' enablegroup 'tag' 'plotgroups' } ...
        {} {'style' 'checkbox'   'string' 'Detach plots' 'value' detachplots 'enable' 'on' 'tag' 'detachtag' } ...
        {} ...
           {'style' 'text'       'string' 'Multiple channels selection' 'fontweight' 'bold' 'tag', 'spec' 'fontsize', 12} ...
        {} {'style' 'radio'      'string' 'Plot channels in scalp array'    'value' radio_scalparray 'tag' 'scalparray'       'userdata' 'radio' 'callback' cb_radio} { } ...
        {} {'style' 'radio'      'string'  'Plot topography at time (ms)' 'value' radio_scalptopo  'tag' 'scalptopotext' 'userdata' 'radio' 'callback' cb_radio} ...
           {'style' 'edit'       'string' num2str(STUDY.etc.erpparams.topotime) 'tag' 'topotime' 'callback' cb_edit } ...
        {} {'style' 'radio'      'string' 'Average selected channels' 'value' radio_averagechan 'tag' 'averagechan' 'userdata' 'radio' 'callback' cb_radio} { } };
    cbline = [0.07 1.1];
    otherline = [ 0.07 0.6 .3];
    chanline  = [ 0.07 0.8 0.3];
    geometry = { 1 otherline otherline otherline 1 1 cbline cbline cbline 1 1 chanline chanline chanline };
    geomvert = [1.2 1 1 1 0.5 1.2 1 1 1 0.5 1.2 1 1 1 ];

    % component plotting
    % ------------------
    if isnan(STUDY.etc.erpparams.topotime)
        geometry(end-4:end) = [];
        geomvert(end-4:end) = [];
        uilist(end-10:end) = [];
    end;

    [out_param userdat tmp res] = inputgui( 'geometry' , geometry, 'uilist', uilist, 'geomvert', geomvert, ...
                                            'title', 'ERP plotting options -- pop_erpparams()');
    if isempty(res), return; end;

    % decode inputs
    % -------------
    %if res.plotgroups & res.plotconditions, warndlg2('Both conditions and group cannot be plotted on the same panel'); return; end;
    if res.plotgroups, res.plotgroups = 'together'; else res.plotgroups = 'apart'; end;
    if res.plotconditions , res.plotconditions  = 'together'; else res.plotconditions  = 'apart'; end;
    if res.detachtag, res.detachtag = 'on'; else res.detachtag = 'off'; end;
    if ~isfield(res, 'topotime'), res.topotime = STUDY.etc.erpparams.topotime;
    else res.topotime  = str2num( res.topotime );
    end;
    res.timerange = str2num( res.timerange );
    res.ylim      = str2num( res.ylim );
    res.filter    = str2num( res.filter );
    if ~isfield(res, 'averagechan'), res.averagechan = STUDY.etc.erpparams.averagechan;
    elseif res.averagechan, res.averagechan = 'on'; else res.averagechan = 'off';
    end;

    % build command call
    % ------------------
    options = {};
    if ~strcmpi( char(res.filter), char(STUDY.etc.erpparams.filter)), options = { options{:} 'filter' res.filter }; end;
    if ~strcmpi( res.plotgroups, STUDY.etc.erpparams.plotgroups), options = { options{:} 'plotgroups' res.plotgroups }; end;
    if ~strcmpi( res.plotconditions , STUDY.etc.erpparams.plotconditions ), options = { options{:} 'plotconditions'  res.plotconditions  }; end;
    if ~strcmpi( res.detachtag, STUDY.etc.erpparams.detachplots), options = { options{:} 'detachplots' res.detachtag}; end;
    if ~isequal(res.ylim       , STUDY.etc.erpparams.ylim),      options = { options{:} 'ylim' res.ylim       }; end;
    if ~isequal(res.timerange  , STUDY.etc.erpparams.timerange)  , options = { options{:} 'timerange' res.timerange }; end;
    if ~isequal(res.averagechan, STUDY.etc.erpparams.averagechan), options = { options{:} 'averagechan' res.averagechan }; end;
    if (all(isnan(res.topotime)) & all(~isnan(STUDY.etc.erpparams.topotime))) | ...
            (all(~isnan(res.topotime)) & all(isnan(STUDY.etc.erpparams.topotime))) | ...
                (all(~isnan(res.topotime)) & ~isequal(res.topotime, STUDY.etc.erpparams.topotime))
        options = { options{:} 'topotime' res.topotime };
    end;
    if ~isempty(options)
        STUDY = pop_erpparams(STUDY, options{:});
        com = sprintf('STUDY = pop_erpparams(STUDY, %s);', vararg2str( options ));
    end;
else
    if strcmpi(varargin{1}, 'default')
        STUDY = default_params(STUDY);
    else
        for index = 1:2:length(varargin)
            if ~isempty(strmatch(varargin{index}, fieldnames(STUDY.etc.erpparams), 'exact'))
                STUDY.etc.erpparams = setfield(STUDY.etc.erpparams, varargin{index}, varargin{index+1});
            end;
        end;
    end;
end;

% scan clusters and channels to remove erpdata info if timerange has changed
% ----------------------------------------------------------
if ~isequal(STUDY.etc.erpparams.timerange, TMPSTUDY.etc.erpparams.timerange)
    rmfields = { 'erpdata' 'erptimes' };
    for iField = 1:length(rmfields)
        if isfield(STUDY.cluster, rmfields{iField})
            STUDY.cluster = rmfield(STUDY.cluster, rmfields{iField});
        end;
        if isfield(STUDY.changrp, rmfields{iField})
            STUDY.changrp = rmfield(STUDY.changrp, rmfields{iField});
        end;
    end;
end;

function STUDY = default_params(STUDY)
    if ~isfield(STUDY.etc, 'erpparams'), STUDY.etc.erpparams = []; end;
    if ~isfield(STUDY.etc.erpparams, 'topotime'),         STUDY.etc.erpparams.topotime = []; end;
    if ~isfield(STUDY.etc.erpparams, 'filter'),           STUDY.etc.erpparams.filter = []; end;
    if ~isfield(STUDY.etc.erpparams, 'timerange'),        STUDY.etc.erpparams.timerange = []; end;
    if ~isfield(STUDY.etc.erpparams, 'ylim'     ),        STUDY.etc.erpparams.ylim      = []; end;
    if ~isfield(STUDY.etc.erpparams, 'plotgroups') ,      STUDY.etc.erpparams.plotgroups = 'apart'; end;
    if ~isfield(STUDY.etc.erpparams, 'plotconditions') ,  STUDY.etc.erpparams.plotconditions  = 'apart'; end;
    if ~isfield(STUDY.etc.erpparams, 'averagechan') ,     STUDY.etc.erpparams.averagechan  = 'off'; end;
    if ~isfield(STUDY.etc.erpparams, 'detachplots') ,     STUDY.etc.erpparams.detachplots  = 'on'; end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function newS = align_shape( TransM , S )

newS = zeros(size(S));

x = S(1:2:end);
y = S(2:2:end);

xy = [x'; y'; ones(1, size(x,1))];
newXY = [TransM ; [0 0 1]]\xy;

newS(1:2:end) = newXY(1,:);
newS(2:2:end) = newXY(2,:);

end

author: tntrung


text:
% takes NimgFirstRegistration mean image and aligns it to itself
% does nonrigid registration across specified blocks
% returns mean image for registration

function ops = nonrigidAlignIterative(data, ops)

% uu = squeeze(sum(sum(data(:,:,:).^2,1),2));
% [~, isort] = sort(uu, 'descend');
% ops.mimg        = data(:,:,isort(50));

% take most correlated frames to compute initial mean image
mimg = pick_reg_init(data);

numBlocks = numel(ops.yBL);
dsold = zeros(size(data,3), 2, numBlocks);
err = zeros(ops.NiterPrealign, 1);
%%
for i = 1:ops.NiterPrealign
    dsnew = zeros(size(data,3), 2, numBlocks,'double');
    Corr = zeros(size(data,3), numBlocks,'double');
    for ib = 1:numBlocks
        % collect ds
        ops.mimg = mimg(ops.yBL{ib},ops.xBL{ib});
        % compute offsets from each frame to mean image and subtract mean
        % offset (centered at offset = 0)
	if ops.kriging
	  [dsnew(:,:,ib), Corr(:,ib)]  = ...
	      regoffKriging(data(ops.yBL{ib},ops.xBL{ib},:), ops, 1);
	else
	  [dsnew(:,:,ib), Corr(:,ib)]  = ...
	      regoffLinear(data(ops.yBL{ib},ops.xBL{ib},:), ops, 1);
	end
    end

    % register frames
    dreg = nonrigidRegFrames(data, ops.xyMask, dsnew);

    % sort by correlation and take mean of most correlated frames
    [~, igood] = sort(mean(Corr,2), 'descend');
    if i<floor(ops.NiterPrealign/2)
        igood = igood(1:100);
    else
        igood = igood(1:round(size(data,3)/2));
    end
    mimg = mean(dreg(:,:,igood),3);

    err(i) = mean(mean(sum((dsold - dsnew).^2,2),3)).^.5;

    dsold = dsnew;
end
ops.mimg = mimg;

ops.AlignNanThresh = median(Corr) - 4*std(Corr);
ops.ErrorInitialAlign = err;
ops.dsprealign = dsnew;

end

author: cortex-lab


text:
function a=recor2(bsol,para)
time=bsol(5)+bsol(6);
kk=trajt(para,bsol);
tt=torque(kk);
ft=ftorque(tt);
ffq=kk(1:2,:);
fq=sum(sum(abs(diff(ffq'))));
pos=forkin(kk(1:2,:));
xx=pos(1,:);yy=pos(2,:);
x=diff(xx);
y=diff(yy);
fdis=sum(sqrt(x.^2+y.^2));
a=[time,fq,fdis,ft];


text:
function [s,it]=display(a,tab,it)

  global display_tree_depth;
  if isempty(display_tree_depth) display_tree_depth=100; end;
  global display_tree_tab_step;
  if isempty(display_tree_tab_step) display_tree_tab_step=3; end;
  global display_tree_show_brackets;
  if isempty(display_tree_show_brackets) display_tree_show_brackets=1; end;

  % ----- tabulation stuff ----------
  print=0;  if nargin==1  tab=0; print=1; it=[]; end;
  tab_step=display_tree_tab_step;
  % ------ display tree ------------
  [s,it]=display(a.algorithm,tab,it);
  if  display_tree_show_brackets
    s=[s ones(1,tab+min(2,tab_step))*32 '{\n' ];
  end
  tab=tab+tab_step;
  i=length(it); v=it{i}; v=[v 1]; it{i}=v; %% add 1 to last row
  for i=1:length(a.child)
     [st it]=display(a.child{i},tab,it);
     if length(it{length(it)})<display_tree_depth
      s =[s ones(1,tab)*32 st];
     end

     j=length(it); v=it{j}; v(length(v))=v(length(v))+1;
     it{j}=v; %% increment last row
  end
  if  display_tree_show_brackets
    if  length(it{length(it)})>=display_tree_depth
     s=[s ones(1,tab)*32 '...\n'];
    end
    s =[s ones(1,tab-tab_step+min(2,tab_step))*32 '}\n' ];
  end
  tab=tab-tab_step;
  %% pop index of last item
  i=length(it); v=it{length(it)}; v=v(1:length(v)-1); it{i}=v;
  if print
   s=s(1:length(s)-2); % remove last carriage return
   it=deal(it(1:length(it)-1)); % remove last index to tree
   disp(sprintf(char(s)));
  end