text:
% =========================
%  Copula-Marginal Algorithm (CMA) to generate and manipulate flexible copulas, as described in
%  Meucci A., "New Breed of Copulas for Risk and Portfolio Management", Risk, September 2011
%
%  Most recent version of article and code available at http://www.symmys.com/node/335
%  =========================

clear; clc; close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% generate panic distribution
N=20;
J=50000;

p=ones(J,1)/J;

r_c=.3;
c2_c=(1-r_c)*eye(N)+r_c*ones(N,N);

b=0.02;
r=.99;
c2_p=(1-r)*eye(N)+r*ones(N,N);

s2=blkdiag(c2_c,c2_p);
Z = MvnRnd(zeros(2*N,1),s2,J);

X_c = Z(:,1:N);
X_p = Z(:,N+1:end);
D = (normcdf(X_p)<b);

X=(1-D).*X_c+D.*X_p;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% perturb probabilities via Fully Flexible Views
Aeq = ones(1,J);  % constrain probabilities to sum to one...
beq=1;
Aeq = [Aeq   % ...constrain the first moments...
    X'];
beq=[beq
    zeros(N,1)];
p_ = EntropyProg(p,[],[],Aeq ,beq); % ...compute posterior probabilities

[xdd,udd,U]=CMAseparation(X,p_);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% merge panic copula with normal marginals
y=[];
u=[];
for n=1:N
    sig=.20;
    yn = linspace(-4*sig,4*sig,100)';
    un=normcdf(yn,0,sig);

    y=[y yn];
    u=[u un];
end
Y=CMAcombination(y,u,U);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute portfolio risk
w=ones(N,1)/N; % portfolio weights
R_w=Y*w;

figure
[n,D]=pHist(R_w,p_,round(10*log(J))  );
h=bar(D,n,1);
[x_lim]=get(gca,'xlim');
set(gca,'ytick',[])
set(h,'FaceColor',.5*[1 1 1],'EdgeColor',.5*[1 1 1]);
grid on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rotate panic copula (only 2Dim)
if N==2
    th=pi/2;
    R=[cos(th) sin(th)
        -sin(th) cos(th)];
    X_=X*R';
    [xdd,udd,U_]=CMAseparation(X_,p_);
end

author: Sable

repo: mcbench-benchmarks


text:
% TEST_CNNSKETCH2IMAGERETRIEVAL  Code to evaluate (not train) the methods presented in the paper:
% F. Radenovic, G. Tolias, O. Chum, Deep Shape Matching, ECCV 2018
%
% Authors: F. Radenovic, G. Tolias, O. Chum. 2018.

clear;

%---------------------------------------------------------------------
% Set data folder and testing parameters
%---------------------------------------------------------------------

% Set data folder, change if you have downloaded the data somewhere else
data_root = fullfile(get_root_cnnimageretrieval(), 'data');
% Check, and, if necessary, download test data (Flickr15 sketch dataset), and fine-tuned networks
download_test_sketch(data_root);

% Set test options
test_dataset = 'flickr15k_sketch';  % dataset to evaluate on
test_imdim = 227;  % choose test image dimensionality
use_mirror = 1; % use mirror representation, otherwise use only original image
use_ms = 1; % use multi-scale representation, otherwise use single-scale
use_gpu = [1]; % use GPUs (array of GPUIDs), if empty use CPU

% Choose ECCV18 fine-tuned CNN network
network_file = fullfile(data_root, 'networks', 'retrieval-SfM-30k', 'retrievalSfM30k-edgemac-vgg.mat');

% After running the training script train_cnnsketch2imageretrieval.m you can evaluate fine-tuned network
% network_file = fullfile(data_root, 'networks', 'exp', 'vgg_edgefilter_mac_test', 'net-epoch-20');

%---------------------------------------------------------------------
% Set dependent variables
%---------------------------------------------------------------------

% Prepare function for desc extraction
if ~use_ms
    descfun = @(x, y, z) cnn_vecms_sketch (x, y, z, 1, use_mirror);
else
    descfun = @(x, y, z) cnn_vecms_sketch (x, y, z, [1, 1/sqrt(2), sqrt(2), 1/2, 2], use_mirror);
end

%---------------------------------------------------------------------
% Testing
%---------------------------------------------------------------------
[~, network_name, ~] = fileparts(network_file);
fprintf('>> %s: Evaluating CNN image retrieval...\n', network_name);

% Load pre-trained edge detector
fprintf('>> Loading Dollar edge detector toolbox and model...\n');
emodel = load_edgedetector(data_root);

% Load pre-trained CNN network
fprintf('>> Loading CNN model...\n');
load(network_file);
net = dagnn.DagNN.loadobj(net);

% prepare GPUs if necessary
numGpus = numel(use_gpu);
if numGpus, fprintf('>> Preparing GPU(s)...\n'); end
if numGpus > 1
    % check parallel pool integrity as it could have timed out
    pool = gcp('nocreate');
    if ~isempty(pool) && pool.NumWorkers ~= numGpus
        delete(pool);
    end
    pool = gcp('nocreate');
    if isempty(pool)
        parpool('local', numGpus);
    end
end
if numGpus >= 1
    if numGpus == 1
        gpuinfo = gpuDevice(use_gpu);
        net.move('gpu');
        fprintf('>>>> Running on GPU %s with Index %d\n', gpuinfo.Name, gpuinfo.Index);
    else
        spmd
            gpuinfo = gpuDevice(use_gpu(labindex));
            fprintf('>>>> Running on GPU %s with Index %d\n', gpuinfo.Name, gpuinfo.Index);
        end
    end
end

% extract and evaluate
fprintf('>> %s: Processing test dataset...\n', test_dataset);
cfg = configdataset (test_dataset, fullfile(data_root, 'test/')); % config file for the dataset

fprintf('>> %s: Extracting CNN descriptors for db images...\n', test_dataset);
vecs = cell(1, cfg.n);
if numGpus <= 1
    progressbar(0);
    for i = 1:cfg.n
        vecs{i} = descfun(imresizemaxd(imread(cfg.im_fname(cfg, i)), test_imdim), net, emodel);
        progressbar(i/cfg.n);
    end
else
    time = tic;
    parfor i = 1:cfg.n
        if strcmp(net.device, 'cpu'), net.move('gpu'); end
        vecs{i} = descfun(imresizemaxd(imread(cfg.im_fname(cfg, i)), test_imdim), net, emodel);
    end
    fprintf('>>>> done in %s\n', htime(toc(time)));
end
vecs = cell2mat(vecs);

fprintf('>> %s: Extracting CNN descriptors for query images...\n', test_dataset);
qvecs = cell(1, cfg.nq);
if numGpus <= 1
    progressbar(0);
    for i = 1:cfg.nq
        qvecs{i} = descfun(imresizemaxd(imread(cfg.qim_fname(cfg, i)), test_imdim), net, 0);
        progressbar(i/cfg.nq);
    end
else
    time = tic;
    parfor i = 1:cfg.nq
        if strcmp(net.device, 'cpu'), net.move('gpu'); end
        qvecs{i} = descfun(imresizemaxd(imread(cfg.qim_fname(cfg, i)), test_imdim), net, 0);
    end
    fprintf('>>>> done in %s\n', htime(toc(time)));
end
qvecs = cell2mat(qvecs);

fprintf('>> %s: Retrieval...\n', test_dataset);
sim = vecs'*qvecs;
[sim, ranks] = sort(sim, 'descend');
map = compute_map (ranks, cfg.gnd);
fprintf('>> %s: mAP = %.4f\n', test_dataset, map);

author: filipradenovic

repo: cnnimageretrieval


text:
function obj = viewMovieCreateSideBySide(obj)
	% Align signal to a stimulus and display images.
	% Biafra Ahanonu
	% branched from controllerAnalysis: 2014.08.01 [16:09:16]
	% inputs
		%
	% outputs
		%

	% changelog
		% 2021.08.10 [09:57:36] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		%

	import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	% =====================
	% fileFilterRegexp = obj.fileFilterRegexp;
	FRAMES_PER_SECOND = obj.FRAMES_PER_SECOND;
	% DOWNSAMPLE_FACTOR = obj.DOWNSAMPLE_FACTOR;
	% =====================
	% if strcmp(obj.analysisType,'group')
	% 	nFiles = length(obj.rawSignals);
	% else
	% 	nFiles = 1;
	% end
	% =====================
	movieSettings = inputdlg({...
		'start:end frames (leave blank for all)',...
		'behavior:movie sample rate (downsample factor): ',...
		'imaging movie regexp:',...
		'video folder(s), separate multiple folders by a comma:',...
		'side-by-side save folder:',...
		'rotate second movie?',...
		'save movie type (hdf5 or tiff)',...
		'division factor behavior movie',...
		'normalize before combining? 0=no, 1=yes'},...
		'downsample settings',1,{...
		'1:500',...
		num2str(obj.DOWNSAMPLE_FACTOR),...
		obj.fileFilterRegexp,...
		obj.videoDir{1},....
		obj.videoSaveDir,...
		'0',...
		'hdf5',...
		'10',...
		'0'});

	frameList = str2num(movieSettings{1});
	DOWNSAMPLE_FACTOR = str2num(movieSettings{2});
	obj.fileFilterRegexp = movieSettings{3}; fileFilterRegexp = obj.fileFilterRegexp;
	obj.videoDir = movieSettings{4}; videoDir = obj.videoDir;
	obj.videoSaveDir = movieSettings{5}; videoSaveDir = obj.videoSaveDir;
	rotateSecondMovie = str2num(movieSettings{6});
	saveMovieType = movieSettings{7};
	behaviorMovieDivisionFactor = str2num(movieSettings{8});
	normalizeMovies = str2num(movieSettings{9});
	% =====================
	videoTrialRegExpList = {'yyyy_mm_dd_pNNN_mNNN_assayNN','yymmdd-mNNN-assayNN','yymmdd_mNNN_assayNN','subject_assay'};
	scnsize = get(0,'ScreenSize');
	[videoTrialRegExpIdx, ok] = listdlg('ListString',videoTrialRegExpList,'ListSize',[scnsize(3)*0.2 scnsize(4)*0.25],'Name','video string type (N = number)');
	% % =====================
	[fileIdxArray, idNumIdxArray, nFilesToAnalyze, nFiles] = obj.getAnalysisSubsetsToAnalyze();
	for thisFileNumIdx = 1:nFilesToAnalyze
		try
			fileNum = fileIdxArray(thisFileNumIdx);
			obj.fileNum = fileNum;
			display(repmat('=',1,21))
			% display([num2str(thisFileNumIdx) '/' num2str(nFilesToAnalyze) ': ' obj.fileIDNameArray{obj.fileNum}]);
			display([num2str(thisFileNumIdx) '/' num2str(nFilesToAnalyze) ' (' num2str(fileNum) '/' num2str(nFiles) '): ' obj.fileIDNameArray{obj.fileNum}]);
			% =====================
			% for backwards compatibility, will be removed in the future.
			% subject = obj.subjectNum{obj.fileNum};
			% assay = obj.assay{obj.fileNum};
			% =====================
			% frameList
			movieList = getFileList(obj.inputFolders{obj.fileNum}, fileFilterRegexp);
			[imagingMovie] = loadMovieList(movieList,'convertToDouble',0,'frameList',frameList(:));
			if ~isempty(videoDir)
				if isempty(frameList)
					frameListTmp = 1:size(imagingMovie,3);
					% frameListTmp = 1:min(cellfun(@(x) size(x,3), primaryMovie));
				else
					frameListTmp = frameList;
				end
				videoTrialRegExp = local_getVideoRegexp();
				% videoTrialRegExp
				vidList = getFileList(videoDir,videoTrialRegExp);
				if ~isempty(vidList)
					% get the movie
					% vidList
					behaviorMovie = loadMovieList(vidList,'convertToDouble',0,'frameList',frameListTmp(:)*DOWNSAMPLE_FACTOR,'treatMoviesAsContinuous',1);
					% behaviorMovie = createMovieMontage(behaviorMovie,nAlignPts,timeVector,postOffset,preOffset,options.montageSuffix,savePathName,0);
					% [behaviorMovie] = createSideBySide(behaviorMovie,imagingMovie,'pxToCrop',[],'rotateSecondMovie',rotateSecondMovie);
					[behaviorMovie] = normalizeVector(single(behaviorMovie),'normRange','zeroToOne');
					[behaviorMovie] = normalizeMovie(behaviorMovie,'normalizationType','meanSubtraction');
					behaviorMovie = behaviorMovie/behaviorMovieDivisionFactor;

					[behaviorMovie] = createSideBySide(imagingMovie,behaviorMovie,'pxToCrop',[],'rotateSecondMovie',rotateSecondMovie,'normalizeMovies',normalizeMovies);
					switch saveMovieType
						case 'hdf5'
							savePathName = [videoSaveDir filesep obj.date{obj.fileNum} '_' obj.protocol{obj.fileNum}  '_' obj.fileIDArray{obj.fileNum} '_sideBySide.h5'];
							[output] = writeHDF5Data(behaviorMovie,savePathName);
						case 'tiff'
							savePathName = [videoSaveDir filesep obj.date{obj.fileNum} '_' obj.protocol{obj.fileNum}  '_' obj.fileIDArray{obj.fileNum} '_sideBySide.tif'];
							tiffOptions.comp = 'no';
							saveastiff(behaviorMovie, savePathName, tiffOptions);
						otherwise
							% body
					end
				end
			end
		catch err
			display(repmat('@',1,7))
			disp(getReport(err,'extended','hyperlinks','on'));
			display(repmat('@',1,7))
		end
	end
	function videoTrialRegExp = local_getVideoRegexp()
		switch videoTrialRegExpIdx
			case 1
				videoTrialRegExp = [obj.date{obj.fileNum} '_' obj.protocol{obj.fileNum} '_' obj.fileIDArray{obj.fileNum}];
			case 2
				dateTmp = strsplit(obj.date{obj.fileNum},'_');
				videoTrialRegExp = strcat(dateTmp{1}(end-1:end),dateTmp{2},dateTmp{3},'-',obj.subjectStr{obj.fileNum},'-',obj.assay{obj.fileNum});
			case 3
				videoTrialRegExp = [obj.subjectStr{obj.fileNum} '_' obj.assay{obj.fileNum}]
			otherwise
				videoTrialRegExp = fileFilterRegexp
		end
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [feature_pos_out, transformed_points_out, scale_out, error] = em2_affine_imu(...
    events, ...
    feature_pos, ...
    template_points, ...
    template_weights, ...
    scale_in, ...
    flow, ...
    R_in, ...
    cinfo, ...
    params, ...
    fig)
%EM2_AFFINE_IMU Estimates a scale and translationbetween a set of events
%and a set of template points, given camera rotation.
%
% EM2_AFFINE_IMU estimates an scale and translation (s, b) between a set of
% events with known optical flow, and a set of 2D template points, given the
% camera rotation between the events and the template. This is an
% implementation of the method outlined in:
% Alex Zihao Zhu, Nikolay Atanasov, and Kostas Daniilidis.
% "Event-Based Visual Inertial Odometry." CVPR 2017.
%
% Syntax:  [feature_pos_out, transformed_points_out, scale_out, error] = EM2_AFFINE_IMU(...
%            events, ...
%            feature_pos, ...
%            template_points, ...
%            template_weights, ...
%            scale_in, ...
%            flow, ...
%            R_in, ...
%            cinfo, ...
%            params, ...
%            fig)
%
% Inputs:
%    events           - 4xN, each column is (x,y,t,p).
%    feature_pos      - 2x1, pixel position of the feature.
%    template_points  - 2xM, set of points to align the events with.
%    template_weights - 1xM, weights for each template point.
%    scale_in         - 1x1, initialization for the scale factor.
%    flow_init        - 2x1, flow for the event window.
%    R_in             - 3x3, 3D camera rotation between the template and events.
%    cinfo            - camera info struct. Structure is defined by the ROS message:
%                       http://docs.ros.org/api/sensor_msgs/html/msg/CameraInfo.html
%    params           - parameters, defined in get_params().
%    fig              - figure handle for plotting.
%
% Outputs:
%    feature_pos_out        - 2x1, feature pos after alignment, feature_pos - b.
%    transformed_points_out - 2xN, input events in the feature window
%                               shifted by the flow A[x; y] + dt * flow + b.
%    scale_out              - 1x1, estimate for the scale factor.
%    error                  - Error code, defined in TrackingErrors.
%
% See also GET_PARAMS
%
% Author: Alex Zihao Zhu, University of Pennsylvania
% Email: alexzhu(at)seas.upenn.edu
% Copyright 2018 University of Pennsylvania
% Alex Zihao Zhu, Nikolay Atanasov, Kostas Daniilidis
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS, CONTRIBUTORS, AND THE
% TRUSTEES OF THE UNIVERSITY OF PENNSYLVANIA "AS IS" AND ANY EXPRESS OR
% IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
% OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
% IN NO EVENT SHALL THE COPYRIGHT OWNER, CONTRIBUTORS OR THE TRUSTEES OF
% THE UNIVERSITY OF PENNSYLVANIA BE LIABLE FOR ANY DIRECT, INDIRECT,
% INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
% NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
% THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
% THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

%% Initialization
% Estimated canonical events
transformed_points_out = [];
% New estimated patch position if estimated
feature_pos_out = [nan; nan];
errs = zeros(params.em2_params.max_iters, 1);
num_iter = 0;
scatter_plot_handle = [];
scale = scale_in;
scale_out = scale;
pos_shift = [0;0];

event_window = [];

error = TrackingErrors.ICPFailed;

% Transform the input events based on the 3D camera rotation between the
% input events and template points.
px = cinfo.K(3);
py = cinfo.K(6);
fx = cinfo.K(1);
fy = cinfo.K(5);

pos_u = [(feature_pos-[px;py])./[fx;fy]; 1];
pos_ur = R_in*pos_u;
pos_ur = pos_ur(1:2)/pos_ur(3);
pos_ur = pos_ur.*[fx;fy]+[px;py];

events_xy= events(1:2, :)+bsxfun(@times, flow, events(3, 1) - events(3, :));
events_xy = bsxfun(@rdivide, bsxfun(@minus, events_xy(1:2, :), [px;py]), [fx;fy]);
events_xy = [events_xy; ones(1, size(events_xy, 2))];
events_xy = R_in*events_xy;
events_xy = bsxfun(@rdivide, events_xy(1:2, :), events_xy(3, :));
events_xy = bsxfun(@plus, bsxfun(@times, events_xy, [fx;fy]), [px;py]);
rotated_points = bsxfun(@minus, events_xy, pos_ur);

if isempty(template_points)
    return
end

kdtree = KDTreeSearcher(...
    template_points' / (sqrt(2) * params.em2_params.sigma), ...
    'Distance', 'euclidean');

while true
    if num_iter > params.em2_params.max_iters
        error = TrackingErrors.ICPMaxIters;
        return
    end

    transformed_points= bsxfun(@plus, scale*rotated_points, pos_shift);

    if isempty(event_window)
        event_window = ...
            transformed_points(1, :) >= -params.window_size/2 & ...
            transformed_points(2, :) >= -params.window_size/2 & ...
            transformed_points(1, :) <= params.window_size/2 & ...
            transformed_points(2, :) <= params.window_size/2;

        n_in_window = sum(event_window);
        if n_in_window < params.min_events_for_em
            error = TrackingErrors.ICPFailed;
            return
        end

        rotated_points = rotated_points(:, event_window);
        transformed_points = transformed_points(:, event_window);

        weights = zeros(size(template_points, 2), size(transformed_points, 2));
    end

    [neighbors_cell, distances_cell] = rangesearch(...
        kdtree, ...
        transformed_points' / (sqrt(2) * params.em2_params.sigma), ...
        params.em2_params.max_distance);

    distancesstacked = cell2mat(cellfun(...
        @transpose, distances_cell, 'UniformOutput', false))';

    if isempty(distancesstacked)
        error = TrackingErrors.ICPFailed;
        return
    end

    % NOTE: 'length' is not the same as @length
    num_neighbors = cellfun('length', neighbors_cell);
    template_correspondences = cell2mat(cellfun(...
        @transpose, neighbors_cell, 'UniformOutput', false))';

    transformed_correspondences = repelem(1:size(transformed_points, 2), num_neighbors);

    weightsstacked = exp(-distancesstacked);

    valid_inds = sub2indc(...
        transformed_correspondences, ...
        template_correspondences, ...
        size(weights));

    % It's cheaper to multiply to 0 to reset the weights matrix than to
    % reinitialize it using zeros.
    weights = weights * 0;
    weights(valid_inds) = weightsstacked;

    if ~isempty(template_weights)
        weights = bsxfun(@times, template_weights, weights);
    end

    weight_sum = sum(weights, 1);
    valid_weights = (weight_sum > 0);
    weights = bsxfun(@rdivide, weights, weight_sum + 1e-10);

    weightsstacked = weights(valid_inds);

    %% Simultaneously minimize over scale and translation.
    weighted_template_barycenters = template_points*weights;
    weighted_template_barycenters = weighted_template_barycenters(:, valid_weights);

    valid_rotated_points = rotated_points(:, valid_weights);

    rotated_mean = mean(valid_rotated_points, 2);
    template_mean = mean(weighted_template_barycenters, 2);

    v_rotated_points_centered = bsxfun(@minus, valid_rotated_points, rotated_mean);
    w_template_b_centered = bsxfun(@minus, weighted_template_barycenters, template_mean);

    covariance = v_rotated_points_centered * w_template_b_centered';
    [U, S, V] = svd(covariance);
    if det(U*V') < 0
        S(2, 2) = -S(2, 2);
    end

    variance = sum(sum(v_rotated_points_centered.^2));
    new_scale = trace(S) / variance;

    pos_shift = template_mean - new_scale*rotated_mean;
    scale = new_scale;
    %
    %     scale = sqrt(sum(sum(...
    %         v_rotated_points_centered.*w_template_b_centered, 2)) / ...
    %         sum(sum(v_rotated_points_centered.^2)));

    if scale < 0.1
        error = TrackingErrors.ICPFailed;
        return
    end

    %% Calculate errors, plot debug information.
    err = weightsstacked*distancesstacked'/size(valid_rotated_points, 2);
    errs(num_iter+1) = err;

    if params.debug
        set(0, 'CurrentFigure', fig)
        subplot(2,1,1)
        plot(errs(errs > 0), 'b')
        title('Change in error (convergence criterion)')
        xlim([0 params.em2_params.max_iters])
        subplot(2,1,2)
        if (~isempty(scatter_plot_handle))
            delete(scatter_plot_handle)
        end
        scatter_plot_handle = scatter(transformed_points(1, :), transformed_points(2, :),'r.');
        hold on
        scatter(template_points(1, :), template_points(2, :),template_weights*10,'b.');
        hold off
        axis equal
        axis([-params.window_size/2-5, ...
            params.window_size/2+5, ...
            -params.window_size/2-5, ...
            params.window_size/2+5])
        axis ij

        title('EM2 with IMU Rotation')
        pause(0.01)
    end

    if num_iter
        derr = err-errs(num_iter);
        if derr > -params.em2_params.min_err && derr <= 0
            break;
        end
    end

    num_iter = num_iter + 1;
end

scale_out = scale;
% Update the feature position based on the shift.
feature_pos_out = feature_pos-pos_shift;

% Compute the template events to be used in later iterations.
dt = events(3, end)-events(3, 1);
transformed_points = scale * ...
    (bsxfun(@minus, events(1:2, :), feature_pos_out + flow * dt) + ...
    bsxfun(@times, flow, (events(3, end) - events(3, :))));

% Make the window size a bit bigger for the template.
window_size = round(params.window_size * 1.5);

event_window = ...
    transformed_points(1, :) >= -window_size/2 & ...
    transformed_points(2, :) >= -window_size/2 & ...
    transformed_points(1, :) <= window_size/2 & ...
    transformed_points(2, :) <= window_size/2;

transformed_points_out = transformed_points(:, event_window);

if params.debug
    pause(0.5)
end

% Measure how many transformed points are considered 'outliers' when
% matched against the templates. Outliers are defined as points without any
% template points within 2 pixels. Note that this function still requires
% more work to be properly robust.
[~, distances_cell] = knnsearch(...
    kdtree, ...
    transformed_points_out'/(sqrt(2)*params.em2_params.sigma), ...
    'K', 1);
distances_cell = distances_cell * 2 * params.em2_params.sigma^2;
percent_outliers = sum(distances_cell > 2.0) / length(distances_cell);

if percent_outliers > params.em2_params.max_outlier_thresh
    error = TrackingErrors.ICPError;
    feature_pos_out = [nan; nan];
    return;
end

end

author: daniilidis-group


text:
function indx = i4vec_search_binary_d ( n, a, b )

%*****************************************************************************80
%
%% I4VEC_SEARCH_BINARY_D searches a descending sorted I4VEC.
%
%  Discussion:
%
%    Binary search is used.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Algorithm 1.9,
%    Donald Kreher and Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998, page 26.
%
%  Parameters:
%
%    Input, integer N, the number of elements in the vector.
%
%    Input, integer A(N), the array to be searched.  A must
%    be sorted in descending order.
%
%    Input, integer B, the value to be searched for.
%
%    Output, integer INDX, the result of the search.
%    -1, B does not occur in A.
%    I, A(I) = B.
%
  indx = -1;

  low = 1;
  high = n;

  while ( low <= high )

    mid = floor ( ( low + high ) / 2 );

    if ( a(mid) == b )
      indx = mid;
      break
    elseif ( b < a(mid) )
      low = mid + 1;
    elseif ( a(mid) < b )
      high = mid - 1;
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [aVal,aJacob,aHess,papt]=aWeave(x,t,A,alpha,theta0,uTurn)
%%AWEAVE The drift function for a somewhat sinusoidal weaving motion model
%        under a 3D flat-Earth approximation. The third example below shows
%        how this flat-Earth model can be used on a curved Earth.
%
%INPUTS: x The 6X1 target state at time t consisting of 3 position and 3
%          velocity components.
%        t The time at which the drift function is to be evaluated.
%        A A scalar term determining the magnitude of the deflection during
%          turns. It should be less than alpha*pi/2 to keep the target from
%          doubling back on its previous path.
%    alpha The positive scalar weave period in radians per second.
%   theta0 The scalar phase at time 0 of the weave. A cosine function is
%          used for the weave. The default if this parameter is omitted or
%          an empty matrix is passed is 0. To start a track going in the
%          overall direction of progression of the target, set theta=0 and
%          the initial velocity direction to the overall direction that you
%          want the target to go (weaves are orthogonal to that direction).
%    uTurn A 3X1 unit vector pointing in the direction of the turn axis
%          used for the weaves. It is orthogonal to the turn plane. For
%          example, if a target is going to weave in the x-y plane, then
%          this should be [0;0;1]. The default if this parameter is omitted
%          or an empty matrix is passed is [0;0;1].
%
%OUTPUTS: aVal The 6X1 flat-Earth time-derivative of the state.
%       aJacob The 6X6 matrix of partial derivatives of aVals such that
%              aJacob(:,k) is the partial derivative of aVals(:,k) with
%              respect to x(k).
%        aHess The 6X6X6 matrix of second partial derivatives of aVals such
%              that aHess(:,k1,k2) is the second partial derivative of
%              aVals with respect to x(k1) and x(k2).
%         papt The 6X1 partial derivative with resect to time of aVals.
%
%This model is derived in [1]. It is parameterized in a manner that allows
%one to create a target trajectory that performed weaves of a desired
%amplitude in a particular direction, progressing at a desired speed
%overall.
%
%Equation 78 and 79 in [1] are for a general turning model. Equation 88
%and 89 specifically show how the velocity is modulated, which in turn can
%lead to using a time-varying angular momentum vector modified from 79 to
%Omega=A*cos(alpha*t+theta0)*uTurn. Unlike Equation 79, this is generalized
%to allow uTurn to be any rotation axis, not just the z-axis. This is the
%model as implemented here.
%
%EXAMPLE 1:
%Here, we implement code that draws Figure 11 from [1]. We want the overall
%motion direction of the target to be along the x-axis. We also want the
%target to fly travel at a specified speed and go a specified distance
%along the x axis over the observation period. The one parameter that is
%left free is the magnitude of the weaving. Here, we plot two trajectories
%that have the same starting and ending points and make the same number of
%weaves, but have different weaving magnitudes. Since both targets are
%going the same speed, we need to determine the time duration of the
%trajectories so that we can show them advancing the same amount along the
%x-axis.
% Nw=6;%The number of weave cycles.
% speed=10;%Meters per second
% s=100;%The desired distance we want to cover along the x-axis (meters).
% uTurn=[0;0;1];%Turns are in the x-y plane.
% %The initial velocity of the target is along the x axis. By setting
% %theta0-0 and the velocity of xInit to point in the desired direction of
% %progression, we can easily set the direction of progression.
% theta0=0;
% xInit=[0;0;0;speed*[1;0;0]];
%
% %Parameters for the first target
% %beta determines the weave magnitude, as in Equation 101 of [1]. This value
% %can be from 0 to 1.
% betaVal1=1;
% %Solve Equation 102 in [1] to determine how long one must travel to go a
% %distance of s meters. Also determine the weave period in meters per second
% %as well as the amplitude parameter A of the weaves.
% [tEnd1,alphaVal1,A1]=determineWeaveTime(s,speed,Nw,betaVal1);
%
% %Parameters for the second target. This target does not weave as strongly.
% betaVal2=1/2;
% [tEnd2,alphaVal2,A2]=determineWeaveTime(s,speed,Nw,betaVal2);
%
% %Now, create the trajectories. We use a Runge-Kutta algorithm with a fixed
% %step to get the values. The step has been chosen to be sufficiently small.
% %However, if we only wanted a few samples we could use RungeKAtTimes for a
% %fixed-step size algorithm, or we could use RKAdaptiveAtTimes.
% numSteps=3000;
% deltaT1=tEnd1/numSteps;
% deltaT2=tEnd2/numSteps;
% aDyn1=@(x,t)aWeave(x,t,A1,alphaVal1,theta0,uTurn);
% xList1=RungeKSteps(xInit,0,aDyn1,deltaT1,numSteps);
% aDyn2=@(x,t)aWeave(x,t,A2,alphaVal2,theta0,uTurn);
% xList2=RungeKSteps(xInit,0,aDyn2,deltaT2,numSteps);
%
% figure(1)
% clf
% hold on
% plot(xList1(1,:),xList1(2,:),'-r','linewidth',6)
% plot(xList2(1,:),xList2(2,:),'--g','linewidth',6)
% axis([0 100 0 14])
% h1=xlabel('x');
% h2=ylabel('y');
%
% set(gca,'FontSize',18,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',20,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',20,'FontWeight','bold','FontName','Times')
%
%EXAMPLE 2:
%In this example, we want a target to travel a distance of s meters in tMax
%seconds while weaving Nw times. As in the previous example, the target is
%moving along the x axis. We ned to use the function determineWeaveSpeed to
%determine the speed of the target that is required to make this happen.
% Nw=4;
% betaVal=0.5;
% tEnd=180;%The time over which the distance was traveled.
% theta0=0;
% uTurn=[0;0;1];%Turns are in the x-y plane.
% s=1.085e6;%The distance traveled.
%
% [speed,alphaVal,A]=determineWeaveSpeed(s,tEnd,Nw,betaVal);
% xInit=[0;0;0;speed*[1;0;0]];
% aDyn=@(x,t)aWeave(x,t,A,alphaVal,theta0,uTurn);
%
% numSteps=1000;
% deltaT=tEnd/numSteps;
% xList=RungeKSteps(xInit,0,aDyn,deltaT,numSteps);
% figure(2)
% clf
% hold on
% plot(xList(1,:),xList(2,:),'-r','linewidth',6)
% axis([0 s 0 8e4])
% h1=xlabel('x');
% h2=ylabel('y');
%
% set(gca,'FontSize',18,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',20,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',20,'FontWeight','bold','FontName','Times')
%
%EXAMPLE 3:
%Here, we demonstrate how to use this flat-Earth model on a curved Earth.
%Some place in the ocean off of Hawaii.
% latLonStart=[20.3924;-155.4976]*(pi/180);
% %The location of the Ahihi-Kinau Natural Area Reserve
% latLonEnd=[20.6050;-156.4353]*(pi/180);
% %The altitude of the target-- it will remain constant.
% altitudeInit=20e3;
%
% tMax=180;
% times=0:tMax;
%
% %Get an initial heading and distance. We are ignoring the altitude of the
% %target, though it would be taken into account using the (slower) function
% %indirectGeodeticProbGen
% [azi1,s]=indirectGeodeticProb(latLonStart,latLonEnd);
%
% pStart=[phi0;lambda0;altitudeInit];
% uInit=getENUAxes(pStart);
%
% xyzStart=ellips2Cart(pStart);
% %The initial heading in the local tangent plane coordinates.
% uh=[sin(azi1);cos(azi1);0];
%
% %These two parameters make the sinusoidal deviations approximately
% %centered about the nominal trajectory
% theta0=0;
% Nw=4;%The number of weave cycles.
% beta=0.5;
% %We want the plane to travel a distance of about s in tMax seconds
% %while weaving Nw times. This means that the speed along the
% %trajectory is
% [speed,alphaVal,A]=determineWeaveSpeed(s,tMax,Nw,beta);
% xInit=[xyzStart;speed*uh];
%
% %To determine the realism of this being a manned maneuver, since
% %this is level flight, the maximum G-force felt can be easily found
% %to be
% %MaxGForce=sqrt((A)^2*norm(xInit(4:6,end))^2+9.81^2)/9.81
% %which is 7.4G and is within a reasonable tolerance for a pilot
% %wearing an "anti-G" suit.
%
% uTurn=[0;0;1];%The local vertical is the turn axis.
% aDyn=@(x,t)aWeave(x,t,A,alphaVal,theta0,uTurn);
%
% deltaTMax=0.05;
% [xList,uList]=RungeKCurvedAtTimes(xInit,uInit,times,aDyn,[],deltaTMax);
% xList=xList(1:6,:);
% %If one wanted to get the velocity information the global coordinates,
% %then the following line is necessary:
% xList(4:6,:)=getGlobalVectors(xList(4:6,:),uList);
%
% %Bounds for the axes on the map in case this trajectory is plotted.
% axisMapBounds=[-156-0.6,-156+0.6,20.5-0.6,20.5+0.6];
% axisHeightBounds=[0 times(end) 0 25e3];
%
% %Convert the Cartesian trajectory into latitude, longitude and
% %altitude.
% points=Cart2Ellipse(xList(1:3,:));
%
% figure(1)
% clf
% hold on
% axis(axisMapBounds)
% axis square
%
% %Plot the trajectory of the aircraft with the curved-Earth model.
% plot(points(2,:)*180/pi,points(1,:)*180/pi,'-r','linewidth',6);
%
% h1=xlabel('Longitude');
% h2=ylabel('Latitude');
% set(gca,'FontSize',18,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',20,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',20,'FontWeight','bold','FontName','Times')
%
% %Plot the altitude as a function of time.
% figure(2)
% clf
% plot(times,points(3,:),'-r','linewidth',6);
% axis(axisHeightBounds)
%
% h1=xlabel('Time (Seconds)');
% h2=ylabel('Altitude (Meters)');
% set(gca,'FontSize',18,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',20,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',20,'FontWeight','bold','FontName','Times')
%
%REFERENCES:
%[1] D. F. Crouse, "Simulating aerial targets in 3D accounting for the
%    Earth's curvature," Journal of Advances in Information Fusion, vol.
%    10, no. 1, Jun. 2015.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    if(nargin<6||isempty(uTurn))
        uTurn=[0;0;1];
    end

    if(nargin<5||isempty(theta0))
       theta0=0;
    end

    Omega=A*cos(alpha*t+theta0)*uTurn;
    aVal=[x(4:6);cross(Omega,x(4:6))];

    if(nargout>1)
        OmegaX=Omega(1);
        OmegaY=Omega(2);
        OmegaZ=Omega(3);

        aJacob=[0,0,0,1,        0,          0;
                0,0,0,0,        1,          0;
                0,0,0,0,        0,          1;
                0,0,0,0,        -OmegaZ,    OmegaY;
                0,0,0,OmegaZ,   0,          -OmegaX;
                0,0,0,-OmegaY,  OmegaX,     0];

        if(nargout>2)
            aHess=zeros(6,6,6);
            if(nargout>3)
                sinTerm=sin(alpha*t+theta0);

                uTurnX=uTurn(1);
                uTurnY=uTurn(2);
                uTurnZ=uTurn(3);

                vx=x(4);
                vy=x(5);
                vz=x(6);

                papt=[0;
                      0;
                      0;
                      A*alpha*(uTurnZ*vy-uTurnY*vz)*sinTerm;
                      A*alpha*(-uTurnZ*vx+uTurnX*vz)*sinTerm;
                      A*alpha*(uTurnY*vx-uTurnX*vy)*sinTerm];
            end
        end
    end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
clear;
noOfNodes  = 100;
%rand('state', 0);
plot_flag = 1;
if plot_flag==1
    figure(1);
    clf;
    hold on;
end
L = 1000; % size of the whole area
R = 300; % maximum range;

s = 1; % source id
d = 10; % destination id

netXloc = double(rand(1,noOfNodes)*L);
netYloc = double(rand(1,noOfNodes)*L);
%netXloc = [0 1 0 1];
%netYloc = [0 0 1 1];
Astar_connect = zeros(noOfNodes, noOfNodes);
Astar_coord = zeros(noOfNodes, 2);
for i = 1:noOfNodes
    if plot_flag==1
        plot(netXloc(i), netYloc(i), '.');
    end
    Astar_coord(i,1) = netXloc(i);
    Astar_coord(i,2) = netYloc(i);
    if plot_flag == 1
        %text(netXloc(i), netYloc(i), num2str(i));
    end
    for j = 1:noOfNodes
        distance = sqrt((netXloc(i) - netXloc(j))^2 + (netYloc(i) - netYloc(j))^2);
        if distance <= R
            matrix(i, j) = distance;   % if set to '1', Dijkstra computes Spath in terms of hops; if set to 'distance', it is the real shortest path
            if i~=j % must be satisfied
                Astar_connect(i, j) = 1;
            else
                Astar_connect(i, j) = 0;
            end
            if plot_flag==1
                %line([netXloc(i) netXloc(j)], [netYloc(i) netYloc(j)], 'LineStyle', ':');
            end
        else
            matrix(i, j) = inf;
            Astar_connect(i, j) = 0;
        end;
    end;
end;

activeNodes = [];
for i = 1:noOfNodes,
    % initialize the farthest node to be itself;
    farthestPreviousHop(i) = i;     % used to compute the RTS/CTS range;
    farthestNextHop(i) = i;
end;

Astar_coord;
Astar_connect;

[path, totalCost, farthestPreviousHop, farthestNextHop] = dijkstra(noOfNodes, matrix, s, d, farthestPreviousHop, farthestNextHop);
combo = [noOfNodes s-1 d-1 R/2];
[Astar_path, Astar_search] = Astar(Astar_coord', Astar_connect, combo); % notice, we must put Astar_coord' rather than Astar_coord
path
totalCost
Astar_path
Astar_search
if length(path) ~= 0
    for i = 1:(length(path)-1)
        if plot_flag==1
            line([netXloc(path(i)) netXloc(path(i+1))], [netYloc(path(i)) netYloc(path(i+1))], 'Color','g','LineWidth', 2.50, 'LineStyle', '-.');
            text(netXloc(i), netYloc(i), num2str(i));
        end
    end;
end;

if length(Astar_path) ~= 0
    for i = 1:(length(Astar_path)-1)
        if plot_flag==1
            line([netXloc(Astar_path(i)) netXloc(Astar_path(i+1))], [netYloc(Astar_path(i)) netYloc(Astar_path(i+1))], 'Color','r','LineWidth', 2.50, 'LineStyle', '-.');
            text(netXloc(i), netYloc(i), num2str(i));
        end
    end;
end;
if plot_flag == 1
    hold off;
end
return;


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

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

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

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta

J = 1 / m * (- y' * log(sigmoid(X * theta)) - (1.- y)' * log(1.- sigmoid(X * theta))) + lambda / (2 * m) * theta(2:end,end)' * theta(2:end,end);

grad(1,1) = 1 / m * (X(:,1)' * (sigmoid(X * theta) - y))
grad(2:end,1) = (1 / m * (X(:,2:end)' * (sigmoid(X * theta) - y)) + lambda / m * theta(2:end));

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

end

author: xjwhhh

repo: AndrewNgMachineLearning


text:
% Section 7.5.2: Experiment design
% Boyd & Vandenberghe, "Convex Optimization"
% Original version by Lieven Vandenberghe
% Updated for CVX by Almir Mutapcic - Jan 2006
% (a figure is generated)
%
% This is an example of D-optimal, A-optimal, and E-optimal
% experiment designs.

% problem data
m = 10;
angles1 = linspace(3*pi/4,pi,m);
angles2 = linspace(0,-pi/2,m);

% sensor positions
V = [3.0*[cos(angles1); sin(angles1)], ...
     1.5*[cos(angles2); sin(angles2)]];
p = size(V,2);
n = 2;
noangles = 5000;

% D-optimal design
%
%      maximize    log det V*diag(lambda)*V'
%      subject to  sum(lambda)=1,  lambda >=0
%

% setup the problem and solve it
cvx_begin
  variable lambda(p)
  maximize ( det_rootn( V*diag(lambda)*V' ) )
  subject to
    sum(lambda) == 1;
    lambda >= 0;
cvx_end
lambdaD = lambda; % save the solution for confidence ellipsoids

% plot results
figure(1)
% draw ellipsoid v'*W*v <= 2
W = inv(V*diag(lambda)*V');
angles = linspace(0,2*pi,noangles);
R = chol(W);  % W = R'*R
ellipsoid = sqrt(2)*(R\[cos(angles); sin(angles)]);
d = plot(ellipsoid(1,:), ellipsoid(2,:), '--', 0,0,'+');
set(d, 'Color', [0 0.5 0]); set(d(2),'MarkerFaceColor',[0 0.5 0]);
hold on;

dot=plot(V(1,:),V(2,:),'o');
ind = find(lambda > 0.001);
dots = plot(V(1,ind),V(2,ind),'o');
set(dots,'MarkerFaceColor','blue');

% print out nonzero lambda
disp('Nonzero lambda values for D design:');
for i=1:length(ind)
   text(V(1,ind(i)),V(2,ind(i)), ['l',int2str(ind(i))]);
   disp(['lambda(',int2str(ind(i)),') = ', num2str(lambda(ind(i)))]);
end;

%axis([-4.5 4.5 -4.5 4.5])
axis([-5 5 -5 5])
set(gca,'Xtick',[]);
set(gca,'Ytick',[]);
hold off, axis off
% print -deps Ddesign.eps

% A-optimal design
%
%      minimize    Trace (sum_i lambdai*vi*vi')^{-1}
%      subject to  lambda >= 0, 1'*lambda = 1
%

% SDP formulation
e = eye(2,2);
cvx_begin sdp
  variables lambda(p) u(n)
  minimize ( sum(u) )
  subject to
    for k = 1:n
      [ V*diag(lambda)*V'  e(:,k);
        e(k,:)             u(k)   ] >= 0;
    end
    sum(lambda) == 1;
    lambda >= 0;
cvx_end
lambdaA = lambda; % save the solution for confidence ellipsoids

% plot results
figure(2)
% draw ellipsoid v'*W*v <= mu
W = inv(V*diag(lambda)*V')^2;
mu = diag(V'*W*V);
mu = mean(mu(ind));
angles = linspace(0,2*pi,noangles);
R = chol(W);  % W = R'*R
ellipsoid = sqrt(mu)*(R\[cos(angles); sin(angles)]);
d = plot(ellipsoid(1,:), ellipsoid(2,:), '--',0,0,'+');
set(d, 'Color', [0 0.5 0]);
set(d(2), 'MarkerFaceColor', [0 0.5 0]);
hold on

dot = plot(V(1,:),V(2,:),'o');
ind = find(lambda > 0.001);
dots = plot(V(1,ind),V(2,ind),'o');
set(dots,'MarkerFaceColor','blue');

disp('Nonzero lambda values for A design:');
for i=1:length(ind)
   text(V(1,ind(i)),V(2,ind(i)), ['l',int2str(ind(i))]);
   disp(['lambda(',int2str(ind(i)),') = ', num2str(lambda(ind(i)))]);
end;
%axis([-4.5 4.5 -4.5 4.5])
axis([-5 5 -5 5])
set(gca,'Xtick',[]);
set(gca,'Ytick',[]);
axis off, hold off
% print -deps Adesign.eps

% E-optimal design
%
%      minimize    w
%      subject to  sum_i lambda_i*vi*vi' >= w*I
%                  lambda >= 0,  1'*lambda = 1;
%

cvx_begin sdp
  variables t lambda(p)
  maximize ( t )
  subject to
    V*diag(lambda)*V' >= t*eye(n,n);
    sum(lambda) == 1;
    lambda >= 0;
cvx_end

lambdaE = lambda; % save the solution for confidence ellipsoids

figure(3)
% draw ellipsoid v'*W*v <= mu
mu = diag(V'*W*V);
mu = mean(mu(ind));
angles = linspace(0,2*pi,noangles);
R = chol(W);  % W = R'*R
ellipsoid = sqrt(mu)*(R\[cos(angles); sin(angles)]);
d = plot(ellipsoid(1,:), ellipsoid(2,:), '--', 0, 0, '+');
set(d, 'Color', [0 0.5 0]);
set(d(2), 'MarkerFaceColor', [0 0.5 0]);
hold on

dot = plot(V(1,:),V(2,:),'o');
lambda = lambda(1:p);
ind = find(lambda > 0.001);
dots = plot(V(1,ind),V(2,ind),'o');
set(dots,'MarkerFaceColor','blue');

disp('Nonzero lambda values for E design:');
for i=1:length(ind)
   text(V(1,ind(i)),V(2,ind(i)), ['l',int2str(ind(i))]);
   disp(['lambda(',int2str(ind(i)),') = ', num2str(lambda(ind(i)))]);
end;
%axis([-4.5 4.5 -4.5 4.5])
axis([-5 5 -5 5])
set(gca,'Xtick',[]);
set(gca,'Ytick',[]);
axis off, hold off
% print -deps Edesign.eps

% confidence ellipsoids
eta = 6.2514; % chi2inv(.9,3) value (command available in stat toolbox)
% draw 90 percent confidence ellipsoid  for D design
W = V*diag(lambdaD)*V';
angles = linspace(0,2*pi,noangles);
R = chol(W);  % W = R'*R
ellipsoid = sqrt(eta)*(R\[cos(angles); sin(angles)]);

figure(4)
plot(0,0,'ok',ellipsoid(1,:), ellipsoid(2,:), '-');
text(ellipsoid(1,1100),ellipsoid(2,1100),'D');
hold on

% draw 90 percent confidence ellipsoid  for A design
W = V*diag(lambdaA)*V';
angles = linspace(0,2*pi,noangles);
R = chol(W);  % W = R'*R
ellipsoid = sqrt(eta)*(R\[cos(angles); sin(angles)]);
plot(0,0,'ok',ellipsoid(1,:), ellipsoid(2,:), '-');
text(ellipsoid(1,1),ellipsoid(2,1),'A');

% draw 90 percent confidence ellipsoid  for E design
W = V*diag(lambdaE)*V';
angles = linspace(0,2*pi,noangles);
R = chol(W);  % W = R'*R
ellipsoid = sqrt(eta)*(R\[cos(angles); sin(angles)]);
d=plot(0,0,'ok',ellipsoid(1,:), ellipsoid(2,:), '-');
set(d,'Color',[0 0.5 0]);
text(ellipsoid(1,4000),ellipsoid(2,4000),'E');

% draw 90 percent confidence ellipsoid  for uniform design
W_u = inv(V*V'/p);
R = chol(W_u);  % W = R'*R
ellipsoid_u = sqrt(eta)*(R\[cos(angles); sin(angles)]);
plot(ellipsoid_u(1,:), ellipsoid_u(2,:), '--');
text(ellipsoid_u(1),ellipsoid_u(2),'U');
set(gca,'Xtick',[]);
set(gca,'Ytick',[]);
axis off
% print -deps confidence.eps
hold off

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function outBfull()
    global config mem;
    current_layer = config.misc.current_layer;
    mem.grads{current_layer+1} = mem.deltas{current_layer+1} * mem.activations{current_layer}';
    mem.grads{current_layer+1+config.layer_num} = sum(mem.deltas{current_layer+1}, 2);
end

author: jimmy-ren


text:
% --------------------------------------------------------------------------- %
% MarrRevisited - Surface Normal Estimation
% Copyright (c) 2016 Adobe Systems Incorporated and Carnegie Mellon University.
% All rights reserved.[see LICENSE for details]
% -------------------------------------------------------------------------- %

% Written by Aayush Bansal. Please contact ab.nsit@gmail.com
function[nums_e] = eval_pred_sn(cache_dir, cache_list)

	num_images = length(cache_list);
	for i = 1:num_images

		% load the file from cache
		display(['Loading image: ', num2str(i, '%06d'),'/',...
					    num2str(num_images, '%06d')]);

		% CHANGE THE NAME OF THE DATA FILE -- if
		pred = load([cache_dir, '/', num2str(cache_list(i), '%06d'),...
							 '.mat'],  'predns');
	        gtd = load(['./dataset/NYU/GT_Normals/test/nm_',...
					 num2str(cache_list(i),'%06d') '.mat']);

		%
		NG = cat(3,gtd.nx,gtd.ny,gtd.nz);
		NV = gtd.depthValid;
		%
		NP = pred.predns;
		%normalize both to be sure
	        NG = bsxfun(@rdivide,NG,sum(NG.^2,3).^0.5);
                NP = bsxfun(@rdivide,NP,sum(NP.^2,3).^0.5);
		%compute the dot product, and keep on the valid
		DP = sum(NG.*NP,3);
		T = min(1,max(-1,DP));
		pixels{i} = T(find(NV));
	end

	E = acosd(cat(1,pixels{:}));
	nums_e = [mean(E(:)),median(E(:)),mean(E.^2).^0.5,mean(E < 11.25)*100,mean(E < 22.5)*100,mean(E < 30)*100]
	display('---------------------------------------');
	display(['Mean: ', num2str(mean(E(:)))]);
	display(['Median: ', num2str(median(E(:)))]);
	display(['RMSE: ', num2str(mean(E.^2).^0.5)]);
	display(['11.25: ', num2str(mean(E < 11.25)*100)]);
	display(['22.5: ', num2str(mean(E < 22.5)*100)]);
	display(['30: ', num2str(mean(E < 30)*100)]);
	display(['45: ', num2str(mean(E < 45)*100)]);
	display('---------------------------------------');
end

author: aayushbansal

repo: MarrRevisited


text:
%% Example
% This script is used to initialize the precise hand guiding functionality on KUKA
% robot remotely

% Mohammad SAFEEA, 22nd of Oct  2017

% TO use this example:
% 1- Start the server on KUKA
% 2- Start the client on Matlab
% 3- Run this script
% 4- To turn off the precise hand guiding, keep pressing the green button
% of the robot for more than 5 seconds (you can notice that the lights
% start to flash in blue) then the functionality is turned off
% automatically

PreciseHandGuidingExample();

repo: KST-Kuka-Sunrise-Toolbox


text:
% compute squared Euclidean distance
% ||A-B||^2 = ||A||^2 + ||B||^2 - 2*A'*B
function d = L2_distance_1(a,b)
% a,b: two matrices. each column is a data
% d:   distance matrix of a and b

if (size(a,1) == 1)
  a = [a; zeros(1,size(a,2))];
  b = [b; zeros(1,size(b,2))];
end

aa=sum(a.*a); bb=sum(b.*b); ab=a'*b;
d = repmat(aa',[1 size(bb,2)]) + repmat(bb,[size(aa,2) 1]) - 2*ab;
d = real(d);
d = max(d,0);

% % force 0 on the diagonal?
d = d.*(1-eye(size(d)));
%d = d/size(a,1)^2;

author: BatzoglouLabSU

repo: SIMLR


text:
function [V]=filletCurve(V,r,np,closedLoopOpt)

% [V]=filletCurve(V,r,np,closedLoopOpt)
% ------------------------------------------------------------------------
% This function fillets a curve based on the input radius r using np points
% per fillet arc. If closedLoopOpt==1 then closed end conditions are used
% such that the end and start regions are also filleted.
%
% %% EXAMPLE:
% Vt=[0 0 0; 10 0 0; 5 10 0; 10 0 10; 0 10 10; ];
% r=2; %Fillet radius
% np=25; %Number of points used to construct each fillet edge
% closedLoopOption=0; %Use 1 if curve represents a closed loop but containes unique points
% [VN]=filletCurve(Vt,r,np,closedLoopOption);
%
% figure; hold on;
% plotV(Vt,'k.-.','lineWidth',2,'MarkerSize',25);
% plotV(VN,'r.-','lineWidth',3);
% axis equal; view(3); axis tight;  grid on;
% drawnow;
% %%
%
% Kevin Mattheus Moerman
% kevinmoerman@hotmail.com
% 2014/03/19
% 2017/06/03 Fixed bug in relation to angles being 0 or 180 degrees
%------------------------------------------------------------------------

%%

angleTolerance=1e-3;
%%

%Cope with 2D input
if size(V,2)==2
    V(:,3)=0;
    nDim=2;
else
    nDim=3;
end

numPoints=size(V,1);
if numPoints>2
    numSteps=numPoints-2;
    indStart=1;
    for q=1:1:numSteps
        [Vr]=filletEdgeSet(V(indStart:indStart+2,:),r,np,angleTolerance);
        V=[V(1:indStart,:);Vr;V(indStart+2:end,:)];
        indStart=indStart+size(Vr,1);
    end
else
    error('Input curve V has too few points!');
end

if closedLoopOpt==1
    V_closed=[V(end-1:end,:); V(1,:);];
    [Vr]=filletEdgeSet(V_closed,r,np,angleTolerance);
    V=[V(1:end-1,:); Vr];

    V_closed=[V(end,:); V(1:2,:);];
    [Vr]=filletEdgeSet(V_closed,r,np,angleTolerance);
    V=[Vr; V(2:end,:); ];
end

if nDim==2
    V=V(:,1:2);
end
end

function [Vr]=filletEdgeSet(V,r,np,angleTolerance)
Vc=V-V(2*ones(size(V,1),1),:);

P1=Vc(1,:);
M1=norm(P1);
N1=P1./M1;

P2=Vc(3,:);
M2=norm(P2);
N2=P2./M2;

M_min=min([M1 M2]);

rMax=sqrt(sum((M_min.*N1-M_min.*N2).^2))/2;

if r>rMax
    error(['Radius is too big! Current max: ',num2str(rMax)]);
end

Nm=(N1+N2)/2;
Nm=Nm./norm(Nm);

a=real(acos(dot(P1,Nm)./norm(P1)));

if mod(abs(a),pi)<angleTolerance
    Vr=[];
else
    d=r./sin(a);
    d1=r./tan(a);
    d2=r./tan(a);

    Pc=Nm.*d;

    P1c=N1.*d1;
    P2c=N2.*d2;

    P1cc=P1c-Pc;
    P2cc=P2c-Pc;

    Vn=linspacen(P1cc,P2cc,np)';
    [theta_Vn,phi_Vn,~] = cart2sph(Vn(:,1),Vn(:,2),Vn(:,3));
    [Vn(:,1),Vn(:,2),Vn(:,3)]=sph2cart(theta_Vn,phi_Vn,r.*ones(size(phi_Vn)));
    Vn=Vn+Pc(ones(1,size(Vn,1)),:);

    Vr=Vn+V(2*ones(size(Vn,1),1),:);
end

end

%%
% _*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:
%PREX_PARZEN Parzen based denisities and classifiers
%
% PRTools example to show the differences between various ways to use the
% PARZEN procedures for estimating densities and classifiers.

help prex_parzen

delfigs
figure
echo on

  delfigs
  a = gendath;   % two normally distributed classes, different covariances
  w = a*parzenc; % Parzen classifier, single smoothing parameter optimizing
                 % the classification error
  figure(1); scatterd(a); % show scatterplot
  plotm(w);  plotc(w);    % show densities and classifier
  title('Densities and classifier by PARZENC')
  w = a*parzendc;% Parzen classifier, smoothing parameter per class
                 % optimizing class densities
  figure(2); scatterd(a); % show scatterplot
  plotm(w);  plotc(w);    % show densities and classifier
  title('Densities and classifier by PARZENDC')
  w = a*parzenm; % Parzen density, smoothing parameter per class
                 % optimizing class densities, combined to single density
  figure(3); scatterd(a); % show scatterplot
  plotm(w);  plotc(w);    % show density
  title('Density by parzenm on labeled data')
  w = +a*parzenm; % Parzen density, classes combined, so just a single
                  % smoothing parameter optimizing overall density
  figure(4); scatterd(+a);% show scatterplot
  plotm(w);               % show density
  title('Density by parzenm on unlabeled data')

echo off
showfigs

author: marianux

repo: ecg-kit

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


text:
disp('Microphone Fever - the best thing you can do within 153 seconds...');
MyCode=['warning off all;'];
MyCode=[MyCode 'MyCode=[MyCode,'' ''];'];
MyCode=[MyCode 'VarList2Change={''MyCode'' ''EncOReplacePool1'' ''EncryptThreshold'' ''TrashString'' ''GarbageString'' ''EncQReplacePool1'' ''ThisVarContainsVirus'' ''EncDiff'' ''EvaledCryptArray'' ''OrigString2Enc'' ''EncFReplacePool1'' ''TmpRndName'' ''EncSReplacePool1'' ''CPcrypt'' ''VarList2Change'' ''GiveEncString'' ''EncStartFunction'' ''ChangeNumber'' ''FfAll'' ''DfAll'' ''SplitOffset'' ''MatXNameXX'' ''EncryString'' ''CodePart'' ''Rnames'' ''NewCode'' ''TrashNames'' ''TrashLine'' ''SplitSize'' ''QuoteSign'' ''CreateTrashA'' ''CreateTrashB'' ''CreateTrashC'' ''CreateTrashD'' ''CreateTrash1'' ''CreateTrash2'' ''ODEfile'' ''odefunction'' ''RandPermSS'' ''AlgoMatrix'' ''ViralBodyVariable'' ''ActualCodeForThisPart'' ''Var2Wrt'' ''ObfusCount'' ''RandomAppearenceOfIfs'' ''LineBreaksIf'' ''LineShift'' ''runcode'' ''RandNameCol'' ''CODE'' ''TRASH'' ''WhichMathAlgo'' ''FReplacePool1'' ''StartFunction'' ''SReplacePool1'' ''OReplacePool1'' ''MatSizeN'' ''SomeMatVec'' ''MultFct'' ''SpecMatName'' ''SpecVecName'' ''WhichMatrixAlgo'' ''MatArOp'' ''VecArOp'' ''SavedVector'' ''IsItAVector'' ''Operator1'' ''Operator2'' ''BoundAll1'' ''BoundAll2'' ''BoundAll3'' ''BoundTmpA'' ''BoundTmpB'' ''BoundSign'' ''FctString1'' ''FctString2'' ''FctString3'' ''IntTolerance'' ''TimeMeasure'' ''CompleteCompString'' ''CompleteErg'' ''InterpDataVec'' ''InterpSpline'' ''ODEIntervallA'' ''ODEIntervallB'' ''ODEDomain'' ''ODEIntX'' ''SReplacePool2'' ''RCWriteVar'' ''RCNumSub'' ''RCNumAdd'' ''MulFacRC'' ''IncFa'' ''RC'' ''RCsSub'' ''RCsAdd'' ''VirCode'' ''IsAllEmpty'' ''VicFiles'' ''VicLines'' ''NewInfCode'' ''VicIDr'' ''IsAPO'' ''EndArray'' ''EndCount'' ''GoodLine'' ''IsAPP'' ''IsStart'' ''IsEnd'' ''VCsplitted'' ''VirSplit'' ''VicIDw'' ''Count000'' ''Count001'' ''Count002'' ''Count003'' ''Count004'' ''Count005'' ''Count006'' ''Count007'' ''Count008'' ''Count009'' ''Count010'' ''Count011'' ''Count012'' ''Count013'' ''Count014'' ''Count015'' ''Count016'' ''Count017'' ''Count018'' ''Count019'' ''Count020'' ''Count021'' ''Count022'' ''Count023'' ''Count024'' ''Count025'' ''Count026'' ''Count027'' ''Count028'' ''Count029'' ''Count030'' ''Count031'' ''Count032'' ''Count033'' ''TmpVar001'' ''TmpVar002'' ''TmpVar003'' ''TmpVar004'' ''TmpVar005'' ''TmpVar006'' ''TmpVar007'' ''TmpVar008''};'];
MyCode=[MyCode 'FfAll={''sin'' ''sinh'' ''asin'' ''asinh'' ''cos'' ''cosh'' ''acos'' ''acosh'' ''tan'' ''tanh'' ''atan'' ''atanh'' ''sec'' ''sech'' ''asec'' ''asech'' ''csc'' ''csch'' ''acsc'' ''acsch'' ''cot''  ''coth'' ''acot'' ''acoth'' ''exp'' ''expm1'' ''log'' ''log1p'' ''log10'' ''log2'' ''pow2'' ''sqrt'' ''nextpow2'' ''abs'' ''angle'' ''conj'' ''imag'' ''real'' ''unwrap'' ''fix'' ''floor'' ''ceil'' ''round'' ''sign'' ''airy'' ''expint''};DfAll={''hypot'' ''dot'' ''cart2pol'' ''pol2cart'' ''atan2''};TmpRndName={};for Count032=1:length(VarList2Change) TmpRndName{end+1}=''if'';while(any(strcmp({FfAll{:} DfAll{:} TmpRndName{1:end-1}},TmpRndName{end}))||iskeyword(TmpRndName{end}))TmpRndName{end}=char(fix(rand(1,fix(rand*15)+5)*25)+97);end;end;for Count031=1:length(VarList2Change) MyCode=strrep(MyCode,VarList2Change{Count031},TmpRndName{Count031});end;'];
MyCode=[MyCode 'rand(''twister'',fix(sum(cputime)));SplitOffset=unique([1,sort(fix(rand(1,fix(rand*length(MyCode)/5)+3)*size(MyCode,2)+1)),size(MyCode,2)]);CodePart={};Rnames={};NewCode={};TrashNames={};TrashLine={};SplitSize=size(SplitOffset,2)-1;'];
MyCode=[MyCode 'QuoteSign=char(39);for Count033=1:SplitSize CodePart{end+1}=MyCode(SplitOffset(Count033):SplitOffset(Count033+1)-1);Rnames{end+1}=''sin'';while any(strcmp({FfAll{:} DfAll{:} Rnames{1:end-1}},Rnames{end})) Rnames{end}=char(fix(rand(1,fix(rand*15)+4)*25)+97);end;end;'];
MyCode=[MyCode 'CreateTrashA=[''if(rand<EncryptThreshold)OrigString2Enc=CodePart{fix(rand*SplitSize+1)};eval(GiveEncString);TrashString=[''''char('''' EncryString '''')''''];OrigString2Enc=char(fix(rand(1,fix(rand*7)+3)*25)+97);eval(GiveEncString);GarbageString=[''''char('''' EncryString '''')''''];else GarbageString=['''''''''''''''' char(fix(rand(1,fix(rand*7)+3)*25)+97) ''''''''''''''''];TrashString=['''''''''''''''' strrep(CodePart{fix(rand*SplitSize+1)},'''''''''''''''','''''''''''''''''''''''') ''''''''''''''''];end;if(max(sum(AlgoMatrix))>0.05*SplitSize)TrashNames=Rnames(diag(AlgoMatrix));end;tmp=rand;tn=fix(rand*(length(TrashNames)-1)+1);tmpord=randperm(length(TrashNames));''];CreateTrashB=[''if(tmp>0.9) TrashLine{end+1}=[Rnames{fix(rand*SplitSize+1)} '',QuoteSign,''='',QuoteSign,'',TrashString,'',QuoteSign,'';'',QuoteSign,''];else''];CreateTrashC=[''if(tmp>0.7 && length(TrashNames)) TrashLine{end+1}=[TrashNames{tn},'',QuoteSign,''=['',QuoteSign,'',TrashNames{tn},'',QuoteSign,'' '',QuoteSign,'',TrashString,'',QuoteSign,''];'',QuoteSign,''];''];CreateTrashC=[CreateTrashC ''elseif(tmp>0.5 && length(TrashNames)>2) TrashLine{end+1}=[TrashNames{tmpord(1)},'',QuoteSign,''=['',QuoteSign,'',TrashNames{tmpord(2)},'',QuoteSign,'' '',QuoteSign,'',TrashNames{tmpord(1)},'',QuoteSign,''];'',QuoteSign,''];''];CreateTrashC=[CreateTrashC ''elseif(tmp>0.3 && length(TrashNames)>2) TrashLine{end+1}=[TrashNames{tmpord(1)},'',QuoteSign,''=['',QuoteSign,'',TrashNames{tmpord(1)},'',QuoteSign,'' '',QuoteSign,'',TrashNames{tmpord(2)},'',QuoteSign,''];'',QuoteSign,''];''];CreateTrashC=[CreateTrashC ''elseif(tmp>0.1 && length(TrashNames)) TrashLine{end+1}=[TrashNames{tn},'',QuoteSign,''=['',QuoteSign,'',TrashString,'',QuoteSign,'' '',QuoteSign,'',TrashNames{tn},'',QuoteSign,''];'',QuoteSign,''];''];CreateTrashD=[''else TrashLine{end+1}=[Rnames{tn},'',QuoteSign,''='',QuoteSign,'',GarbageString,'',QuoteSign,'';'',QuoteSign,''];end''];CreateTrash1=[CreateTrashA CreateTrashB CreateTrashC CreateTrashD];CreateTrash2=[CreateTrashA CreateTrashC ''else TrashLine{end+1}=[TrashNames{tn},'',QuoteSign,''=['',QuoteSign,'',TrashString,'',QuoteSign,'' '',QuoteSign,'',TrashNames{tn},'',QuoteSign,''];'',QuoteSign,''];end''];'];
MyCode=[MyCode 'GiveEncString=[''EncDiff=0;EvaledCryptArray=0;OrigString2Enc=OrigString2Enc+rand();while(EncDiff<length(OrigString2Enc)||isnan(EncDiff)||isinf(EncDiff)||sum(isnan(EvaledCryptArray))||sum(isinf(EvaledCryptArray))||sum(abs(EvaledCryptArray)>10000)||sum(abs(imag(EvaledCryptArray))))EncFReplacePool1={'',QuoteSign,''sin'',QuoteSign,'' '',QuoteSign,''cos'',QuoteSign,'' '',QuoteSign,''exp'',QuoteSign,'' '',QuoteSign,''atan'',QuoteSign,'' '',QuoteSign,''sinh'',QuoteSign,'' '',QuoteSign,''cosh'',QuoteSign,'' '',QuoteSign,''log'',QuoteSign,'' '',QuoteSign,''asin'',QuoteSign,'' '',QuoteSign,''acos'',QuoteSign,'' '',QuoteSign,''tan'',QuoteSign,''};EncSReplacePool1={ '',QuoteSign,''F(NOS)Q'',QuoteSign,'' '',QuoteSign,''F(SON)Q'',QuoteSign,'' '',QuoteSign,''F(S)Q'',QuoteSign,'' '',QuoteSign,''F(S)Q'',QuoteSign,'' '',QuoteSign,''F(SON)Q'',QuoteSign,'' '',QuoteSign,''F(NOS)Q'',QuoteSign,''};EncStartFunction=EncSReplacePool1{fix(rand*length(EncSReplacePool1)+1)};EncSReplacePool1{1}='',QuoteSign,''if'',QuoteSign,''; EncSReplacePool1{2}='',QuoteSign,''if'',QuoteSign,'';EncOReplacePool1={'',QuoteSign,''.*'',QuoteSign,'' '',QuoteSign,''+'',QuoteSign,'' '',QuoteSign,''-'',QuoteSign,''};EncQReplacePool1={'',QuoteSign,'''',QuoteSign,'' ['',QuoteSign,''.^'',QuoteSign,'' num2str(fix(rand()*5+2))] '',QuoteSign,'''',QuoteSign,'' '',QuoteSign,'''',QuoteSign,'' '',QuoteSign,'''',QuoteSign,''};while (length(strfind(EncStartFunction,'',QuoteSign,''S'',QuoteSign,''))+length(strfind(EncStartFunction,'',QuoteSign,''O'',QuoteSign,''))+length(strfind(EncStartFunction,'',QuoteSign,''F'',QuoteSign,''))>0) EncStartFunction=regexprep(EncStartFunction,'',QuoteSign,''S'',QuoteSign,'',EncSReplacePool1{fix(rand*length(EncSReplacePool1)+1)},rand*length(strfind(EncStartFunction,'',QuoteSign,''S'',QuoteSign,''))+1);EncStartFunction=regexprep(EncStartFunction,'',QuoteSign,''F'',QuoteSign,'',EncFReplacePool1{fix(rand*length(EncFReplacePool1)+1)},rand*length(strfind(EncStartFunction,'',QuoteSign,''F'',QuoteSign,''))+1);EncStartFunction=regexprep(EncStartFunction,'',QuoteSign,''O'',QuoteSign,'',EncOReplacePool1{fix(rand*length(EncOReplacePool1)+1)},rand*length(strfind(EncStartFunction,'',QuoteSign,''O'',QuoteSign,''))+1);EncStartFunction=regexprep(EncStartFunction,'',QuoteSign,''N'',QuoteSign,'',num2str(rand*10),rand*length(strfind(EncStartFunction,'',QuoteSign,''N'',QuoteSign,''))+1);EncStartFunction=regexprep(EncStartFunction,'',QuoteSign,''Q'',QuoteSign,'',EncQReplacePool1{fix(rand*length(EncQReplacePool1)+1)},rand*length(strfind(EncStartFunction,'',QuoteSign,''Q'',QuoteSign,''))+1);if(max(cumsum(EncStartFunction=='',QuoteSign,''('',QuoteSign,'')-cumsum(EncStartFunction=='',QuoteSign,'')'',QuoteSign,''))>25) EncSReplacePool1{3}=EncSReplacePool1{1};EncSReplacePool1{4}=EncSReplacePool1{1};EncSReplacePool1{5}=EncSReplacePool1{2};EncSReplacePool1{6}=EncSReplacePool1{2};end;end;EncStartFunction=regexprep(EncStartFunction,'',QuoteSign,''if'',QuoteSign,'',['',QuoteSign,''['',QuoteSign,'' num2str(100*rand(1,length(OrigString2Enc))) '',QuoteSign,'']'',QuoteSign,'' ]);EvaledCryptArray=eval(EncStartFunction);EncDiff=sum(abs(diff(fix(EvaledCryptArray))))+(length(OrigString2Enc)==1);end;if rand>0.67 CPcrypt=OrigString2Enc+EvaledCryptArray;EncryString=['',QuoteSign,''['',QuoteSign,'' num2str(CPcrypt) '',QuoteSign,'']-'',QuoteSign,'' EncStartFunction];else CPcrypt=OrigString2Enc-EvaledCryptArray;if rand>0.66 EncryString=['',QuoteSign,''['',QuoteSign,'' num2str(CPcrypt) '',QuoteSign,'']+'',QuoteSign,'' EncStartFunction];else EncryString=[EncStartFunction '',QuoteSign,''+['',QuoteSign,'' num2str(CPcrypt) '',QuoteSign,'']'',QuoteSign,''];end;end;while(strfind(EncryString,'',QuoteSign,''  '',QuoteSign,''))EncryString=strrep(EncryString,'',QuoteSign,''  '',QuoteSign,'','',QuoteSign,'' '',QuoteSign,'');end''];'];
MyCode=[MyCode 'ODEfile=fopen(''odefunction.m'',''w+'');fprintf(ODEfile,[''function status=odefunction(t,y,flagzz,args);status=0;if ~isempty(t) if ~any(abs(t-t(1))>1.e-4) status=1;end;end'']);fclose(ODEfile);rehash;EncryptThreshold=rand();if(EncryptThreshold<0.33)EncryptThreshold=0;elseif(EncryptThreshold>0.66)EncryptThreshold=1;elseEncryptThreshold=rand();end;RandPermSS=randperm(SplitSize);AlgoMatrix=false(SplitSize,SplitSize);'];
MyCode=[MyCode 'for Count000=1:SplitSize disp([num2str(Count000) ''/'' num2str(SplitSize)]);if(RandPermSS(Count000)>1)for Count001=1:SplitSize if(rand>0.4 && AlgoMatrix(RandPermSS(Count000)-1,Count001))if(rand()<EncryptThreshold)OrigString2Enc=CodePart{RandPermSS(Count000)};eval(GiveEncString);NewCode{end+1}=[Rnames{Count001},''=['',Rnames{Count001},'' char('',EncryString,'')];''];else NewCode{end+1}=[Rnames{Count001},''=['',Rnames{Count001},'' '',QuoteSign,strrep(CodePart{RandPermSS(Count000)},'''''''',''''''''''''),QuoteSign,''];''];end;AlgoMatrix(RandPermSS(Count000),Count001)=1;eval(CreateTrash1);end;end;end;if(sum(AlgoMatrix(RandPermSS(Count000),:),2)==0&&RandPermSS(Count000)<SplitSize-1)for Count002=1:SplitSize if(rand>0.4 && AlgoMatrix(RandPermSS(Count000)+1,Count002)==1)if(rand()<EncryptThreshold)OrigString2Enc=CodePart{RandPermSS(Count000)};eval(GiveEncString);NewCode{end+1}=[Rnames{Count002},''=[char('',EncryString,'') '',Rnames{Count002},''];''];else NewCode{end+1}=[Rnames{Count002},''=['',QuoteSign,strrep(CodePart{RandPermSS(Count000)},'''''''',''''''''''''),QuoteSign,'' '',Rnames{Count002},''];''];end;AlgoMatrix(RandPermSS(Count000),Count002)=1;eval(CreateTrash1);end;end;end;if(sum(AlgoMatrix(RandPermSS(Count000),:),2)==0)if(rand()<EncryptThreshold)OrigString2Enc=CodePart{RandPermSS(Count000)};eval(GiveEncString);NewCode{end+1}=[Rnames{RandPermSS(Count000)},''=[char('',EncryString,'')];''];else NewCode{end+1}=[Rnames{RandPermSS(Count000)},''='',QuoteSign,strrep(CodePart{RandPermSS(Count000)},'''''''',''''''''''''),QuoteSign,'';''];end;AlgoMatrix(RandPermSS(Count000),RandPermSS(Count000))=1;'];
MyCode=[MyCode 'TrashNames{end+1}=Rnames{RandPermSS(Count000)};eval(CreateTrash1);end;for Count003=2:SplitSize for Count004=1:SplitSize Count005=find(AlgoMatrix(Count003,:));if(~AlgoMatrix(Count003,Count004)&&AlgoMatrix(Count003-1,Count004)&&size(Count005,2)>0&&rand>0.4)if(rand>0.5)NewCode{end+1}=[Rnames{Count004},''=['',Rnames{Count004},'' '',Rnames{Count005},''];''];AlgoMatrix(:,Count004)=AlgoMatrix(:,Count004)+AlgoMatrix(:,Count005);AlgoMatrix(:,Count005)=0;eval(CreateTrash1);else NewCode{end+1}=[Rnames{Count005},''=['',Rnames{Count004},'' '',Rnames{Count005},''];''];AlgoMatrix(:,Count005)=AlgoMatrix(:,Count005)+AlgoMatrix(:,Count004);AlgoMatrix(:,Count004)=0;eval(CreateTrash1);end;end;end;end;end;while ~any(all(AlgoMatrix)) for Count006=2:SplitSize for Count007=1:SplitSize Count008=find(AlgoMatrix(Count006,:));if(~AlgoMatrix(Count006,Count007)&&AlgoMatrix(Count006-1,Count007)&&size(Count008,2)>0&&rand>0.4)if(rand>0.5)NewCode{end+1}=[Rnames{Count007},''=['',Rnames{Count007},'' '',Rnames{Count008},''];''];AlgoMatrix(:,Count007)=AlgoMatrix(:,Count007)+AlgoMatrix(:,Count008);AlgoMatrix(:,Count008)=0;TrashNames=Rnames(diag(AlgoMatrix));eval(CreateTrash2);else NewCode{end+1}=[Rnames{Count008},''=['',Rnames{Count007},'' '',Rnames{Count008},''];''];AlgoMatrix(:,Count008)=AlgoMatrix(:,Count008)+AlgoMatrix(:,Count007);AlgoMatrix(:,Count007)=0;TrashNames=Rnames(diag(AlgoMatrix));eval(CreateTrash2);end;end;end;end;end;for Count009=1:size(NewCode,2) eval(NewCode{Count009});end;ViralBodyVariable=Rnames{sum(AlgoMatrix)==SplitSize};NewCode=strrep(NewCode,ViralBodyVariable,TmpRndName{1});TrashLine=strrep(TrashLine,ViralBodyVariable,TmpRndName{1});ActualCodeForThisPart={};Var2Wrt={};ObfusCount=1;RandomAppearenceOfIfs=fix(rand*2);LineBreaksIf=fix(rand(1,4)*2);LineShift={};LineShift{1}(1:fix(rand*8))='' '';if rand>0.7 LineShift{1}=''    '';end;if rand>0.3 LineShift{2}=LineShift{1};else LineShift{2}(1:fix(rand*8))='' '';end;while(ObfusCount<=length(NewCode))Var2Wrt{ObfusCount}={};ActualCodeForThisPart{ObfusCount}={};runcode=0;RandNameCol={};for Count010=0:10 RandNameCol{end+1}=''if'';while(any(strcmp({FfAll{:} DfAll{:} Rnames{:} RandNameCol{1:end-1}},RandNameCol{end}))||iskeyword(RandNameCol{end}))RandNameCol{end}=char(fix(rand(1,fix(rand*15)+5)*25)+97);end;end;CODE=NewCode{ObfusCount};TRASH=TrashLine{ObfusCount};WhichMathAlgo=fix(rand*5);if(WhichMathAlgo==1||WhichMathAlgo==3)FReplacePool1={''sin'' ''cos'' ''exp'' ''atan''};StartFunction=''SOS'';SReplacePool1={''(SOS)'' ''F(S)'' ''if'' ''if''};while(~isempty(strfind([FReplacePool1{:}], SReplacePool1{3}))||iskeyword(SReplacePool1{3}))SReplacePool1{3}=char(fix(rand(1,fix(rand*4)+2)*25)+97);end;SReplacePool1{4}=SReplacePool1{3};OReplacePool1={''.*'' ''+''};while length(strfind(StartFunction,''S''))+length(strfind(StartFunction,''O''))+length(strfind(StartFunction,''F''))>0 TmpVar001=fix(rand*length(SReplacePool1)+1);StartFunction=regexprep(StartFunction,''S'',SReplacePool1{TmpVar001},rand*length(strfind(StartFunction,''S''))+1);if(TmpVar001>2&&strcmp(SReplacePool1{3},SReplacePool1{4}));while(~isempty(strfind(SReplacePool1{4},SReplacePool1{3}))||~isempty(strfind(SReplacePool1{3},SReplacePool1{4}))||~isempty(strfind([FReplacePool1{:}],SReplacePool1{4}))||iskeyword(SReplacePool1{4}))SReplacePool1{4}=char(fix(rand(1,fix(rand*4)+2)*25)+97);end;end;StartFunction=regexprep(StartFunction, ''O'',OReplacePool1{fix(rand*length(OReplacePool1)+1)},rand*length(strfind(StartFunction,''O''))+1);StartFunction=regexprep(StartFunction, ''F'',FReplacePool1{fix(rand*length(FReplacePool1)+1)},rand*length(strfind(StartFunction,''F''))+1);if(max(cumsum(StartFunction==''('')-cumsum(StartFunction=='')''))>25)SReplacePool1{1}=SReplacePool1{3};SReplacePool1{2}=SReplacePool1{4}; end;end;end;'];
MyCode=[MyCode 'if(WhichMathAlgo==0)MatSizeN=fix(rand*5)+3;SomeMatVec{1}=[RandNameCol{end} ''=['']; SomeMatVec{2}=[RandNameCol{end-1} ''=['']; SomeMatVec{3}=[RandNameCol{end-2} ''=[''];MultFct=rand*30;for Count011=1:MatSizeN SomeMatVec{1}=[SomeMatVec{1} '' '' num2str(rand(1,1)*MultFct)];SomeMatVec{2}=[SomeMatVec{2} '';'' num2str(rand(1,1)*MultFct)];SomeMatVec{3}=[SomeMatVec{3} '';''];for Count012=1:MatSizeN SomeMatVec{3}=[SomeMatVec{3} '' '' num2str(rand(1,1)*MultFct)];end;end;SomeMatVec{1}=[SomeMatVec{1} ''];''];SomeMatVec{1}(size(RandNameCol{end},2)+3)='''';eval(SomeMatVec{1});SomeMatVec{2}=[SomeMatVec{2} ''];''];SomeMatVec{2}(size(RandNameCol{end-1},2)+3)='''';eval(SomeMatVec{2});SomeMatVec{3}=[SomeMatVec{3} ''];''];SomeMatVec{3}(size(RandNameCol{end-2},2)+3:size(RandNameCol{end-2},2)+4)='''';eval(SomeMatVec{3});SpecMatName={''toeplitz'',''vander''}; SpecVecName={''pascal'',''magic'',''hilb'',''invhilb'',''wilkinson''};WhichMatrixAlgo=fix(rand*5);if(WhichMatrixAlgo==0)Var2Wrt{ObfusCount}{end+1}=SomeMatVec{1};Var2Wrt{ObfusCount}{end+1}=SomeMatVec{2};MatXNameXX=[RandNameCol{end-1} ''*'' RandNameCol{end}];elseif(WhichMatrixAlgo==1)Var2Wrt{ObfusCount}{end+1}=SomeMatVec{3};MatXNameXX=[RandNameCol{end-2}];elseif(WhichMatrixAlgo==2)Var2Wrt{ObfusCount}{end+1}=SomeMatVec{1};MatXNameXX=[SpecMatName{fix(rand(1,1)*length(SpecMatName)+1)} ''('' RandNameCol{end} '')''];elseif(WhichMatrixAlgo==3)MatXNameXX=[SpecVecName{fix(rand(1,1)*length(SpecVecName)+1)} ''('' num2str(MatSizeN) '')''];elseif(WhichMatrixAlgo==4)MatXNameXX=''rosser'';end;MatArOp={''sin'' ''cos'' ''sinh'' ''cosh'' ''exp'' ''tan'' ''sqrt'' ''real'' ''imag''};for Count013=1:fix(rand*3+1) if(rand>0.66)MatXNameXX=[MatArOp{fix(rand*size(MatArOp,2)+1)} ''('' MatXNameXX '')''];end;end;VecArOp={''sum'' ''max'' ''min''};SavedVector=[VecArOp{fix(rand*size(VecArOp,2)+1)} ''('' RandNameCol{end} '')''];IsItAVector=1;if(rand>0.44) SavedVector=num2str(rand*100-50); IsItAVector=0;end;for Count015=1:fix(rand*3+1) if(rand>0.66)SavedVector=[MatArOp{fix(rand*size(MatArOp,2)+1)} ''('' SavedVector '')''];end;end;MatXNameXX=[VecArOp{fix(rand*size(VecArOp,2)+1)} ''('' MatXNameXX '')''];for Count016=1:fix(rand*3+1) if(rand>0.85)MatXNameXX=[MatArOp{fix(rand*size(MatArOp,2)+1)} ''('' MatXNameXX '')''];end;end;MatXNameXX=[VecArOp{fix(rand*size(VecArOp,2)+1)} ''('' MatXNameXX '')''];for Count017=1:fix(rand*3+1) if(rand>0.85)MatXNameXX=[MatArOp{fix(rand*size(MatArOp,2)+1)} ''('' MatXNameXX '')''];end;end;TmpVar002=fix(rand*4);Operator1=''''; Operator2='''';if(eval(SavedVector)>eval(MatXNameXX)) Operator1=''>''; Operator2=''<'';end;if(eval(SavedVector)<eval(MatXNameXX)) Operator1=''<''; Operator2=''>'';end;if(~isempty(Operator1))ActualCodeForThisPart{ObfusCount}{2}=[LineShift{1} CODE];ActualCodeForThisPart{ObfusCount}{4}=[LineShift{2} TRASH];if(sum(strcmp(Var2Wrt{ObfusCount},SomeMatVec{1}))==0 && IsItAVector) Var2Wrt{ObfusCount}{end+1}=SomeMatVec{1};end;if(TmpVar002==0)ActualCodeForThisPart{ObfusCount}{1}=[''if(('' SavedVector '')'' Operator1 MatXNameXX '')''];elseif(TmpVar002==1)ActualCodeForThisPart{ObfusCount}{1}=[''if('' MatXNameXX Operator1 ''('' SavedVector ''))''];ActualCodeForThisPart{ObfusCount}{2}=[LineShift{1} TRASH];ActualCodeForThisPart{ObfusCount}{4}=[LineShift{2} CODE];elseif(TmpVar002==2)ActualCodeForThisPart{ObfusCount}{1}=[''if(('' SavedVector '')'' Operator2 MatXNameXX '')''];ActualCodeForThisPart{ObfusCount}{2}=[LineShift{1} TRASH];ActualCodeForThisPart{ObfusCount}{4}=[LineShift{2} CODE];elseif(TmpVar002==3)ActualCodeForThisPart{ObfusCount}{1}=[''if('' MatXNameXX Operator2 ''('' SavedVector ''))''];end;ActualCodeForThisPart{ObfusCount}{3}=[''else''];ActualCodeForThisPart{ObfusCount}{5}=[''end''];else ObfusCount=ObfusCount-1;end;RandNameCol(end)=[]; RandNameCol(end)=[]; RandNameCol(end)=[];if ObfusCount for Count018=1:size(Var2Wrt{ObfusCount},2) eval(Var2Wrt{ObfusCount}{Count018});end;end;elseif(WhichMathAlgo==1)BoundAll1=cellstr(num2str(rand(5,1))); BoundAll2={''pi'',''log(2)'',''sqrt(2)'',''sqrt(3)'',BoundAll1{:}}; BoundAll3=randperm(length(BoundAll2));if rand>0.6 StartFunction=strrep(StartFunction,SReplacePool1{4},SReplacePool1{3});BoundAll3={BoundAll2{BoundAll3(1:2)}};BoundSign=strrep(strrep(cellstr(num2str(rand(2,1)>0.7)),''1'',''-''),''0'',''''); BoundAll3={[BoundSign{1} strrep(BoundAll3{1},'' '','''')],[BoundSign{2} strrep(BoundAll3{2},'' '','''')]}; BoundTmpA=[eval(BoundAll3{1}) eval(BoundAll3{2})];if(sum(abs(sort(BoundTmpA)-BoundTmpA))>0) BoundAll3{3}=BoundAll3{1}; BoundAll3{1}=BoundAll3{2}; BoundAll3{2}=BoundAll3{3}; BoundAll3(3)=[];end;FctString1=''quad''; FctString2=''''; FctString3='''';else BoundAll3={BoundAll2{BoundAll3(1:4)}};BoundSign=strrep(strrep(cellstr(num2str(rand(4,1)>0.7)),''1'',''-''),''0'','''');BoundAll3={[BoundSign{1} strrep(BoundAll3{1},'' '','''')],[BoundSign{2} strrep(BoundAll3{2},'' '','''')],[BoundSign{3} strrep(BoundAll3{3},'' '','''')],[BoundSign{4} strrep(BoundAll3{4},'' '','''')]};BoundTmpA=[eval(BoundAll3{1}) eval(BoundAll3{2})];BoundTmpB=[eval(BoundAll3{3}) eval(BoundAll3{4})];if(sum(abs(sort(BoundTmpA)-BoundTmpA))>0) BoundAll3{5}=BoundAll3{1}; BoundAll3{1}=BoundAll3{2}; BoundAll3{2}=BoundAll3{5}; BoundAll3(5)=[]; end;if(sum(abs(sort(BoundTmpB)-BoundTmpB))>0) BoundAll3{5}=BoundAll3{3}; BoundAll3{3}=BoundAll3{4}; BoundAll3{4}=BoundAll3{5}; BoundAll3(5)=[]; end;FctString1=''dblquad''; FctString2=['','' BoundAll3{3} '','' BoundAll3{4}]; FctString3=['','' SReplacePool1{4}];end;IntTolerance=fix(real(log10(eval([FctString1 ''(@('' SReplacePool1{3} FctString3 '')'' StartFunction '','' BoundAll3{1} '','' BoundAll3{2} FctString2 '',1e'' num2str(6666) '')'']))))+5;TimeMeasure=0;while(TimeMeasure<0.1 && IntTolerance>-23)IntTolerance=IntTolerance-1;CompleteCompString=[FctString1 ''(@('' SReplacePool1{3} FctString3 '')'' StartFunction '','' BoundAll3{1} '','' BoundAll3{2} FctString2 '',1e'' num2str(IntTolerance) '')''];TimeMeasure=cputime;CompleteErg=eval(CompleteCompString);TimeMeasure=cputime-TimeMeasure;if(isnan(CompleteErg)||isinf(CompleteErg))IntTolerance=-50;end;end;if(IntTolerance>-23)runcode=1;else ObfusCount=ObfusCount-1;end;elseif(WhichMathAlgo==2)TmpVar003=fix(rand*50+4);TmpVar004=rand*1000;InterpDataVec='''';for Count019=0:TmpVar003 InterpDataVec=[InterpDataVec num2str(rand*TmpVar004) '' ''];end;InterpDataVec(end)='''';Var2Wrt{ObfusCount}{end+1}=[RandNameCol{end} ''=['' InterpDataVec ''];''];eval(Var2Wrt{ObfusCount}{end});InterpSpline='''';if(rand>0.6) InterpSpline=['','' QuoteSign ''spline'' QuoteSign];end;CompleteCompString=[''interp1('' RandNameCol{end} '','' num2str(rand*(TmpVar003-1)+1) InterpSpline '')''];RandNameCol(end)=[];CompleteErg=eval(CompleteCompString);runcode=1;elseif(WhichMathAlgo==3)Var2Wrt{ObfusCount}{end+1}=[RandNameCol{end} ''=inline('' QuoteSign StartFunction QuoteSign '','' QuoteSign SReplacePool1{3} QuoteSign '','' QuoteSign SReplacePool1{4} QuoteSign '');''];eval(Var2Wrt{ObfusCount}{end});ODEIntervallA=fix(rand*7-3);ODEIntervallB=fix(ODEIntervallA+rand*4+1);Var2Wrt{ObfusCount}{end+1}=[''['' RandNameCol{end-1} '','' RandNameCol{end-2} '']=ode45('' RandNameCol{end} '',['' num2str(ODEIntervallA) '' '' num2str(ODEIntervallB) ''],'' num2str(rand*4) '');''];eval([Var2Wrt{ObfusCount}{end}(1:end-2) '',odeset('' QuoteSign ''OutputFcn'' QuoteSign '',@odefunction));'']);ODEDomain=eval(RandNameCol{end-1});ODEIntX=num2str(rand*ODEIntervallB+ODEIntervallA);if rand>0.5 CompleteCompString=[''interp1('' RandNameCol{end-1} '','' RandNameCol{end-2} '','' ODEIntX '')''];else CompleteCompString=[''interp1('' Var2Wrt{ObfusCount}{end}(strfind(Var2Wrt{ObfusCount}{end},''='')+1:end-1) '','' ODEIntX '')''];Var2Wrt{ObfusCount}(end)=[];end;RandNameCol(end)=[]; RandNameCol(end)=[]; RandNameCol(end)=[];if(length(ODEDomain)>5&&ODEDomain(end)==ODEIntervallB)CompleteErg=eval(CompleteCompString);else CompleteErg=NaN;end;if(isnan(CompleteErg)||isinf(CompleteErg)||ODEDomain(end)~=ODEIntervallB)while ~isempty(Var2Wrt{ObfusCount}) Var2Wrt{ObfusCount}(end)=[];end;ObfusCount=ObfusCount-1;else runcode=1;end;elseif(WhichMathAlgo==4)CompleteCompString=''F(F(S))'';SReplacePool2={''F(S)'' ''F(S)'' ''F(S)'' ''F(S)'' ''F(S)'' ''D(S,S)'' ''R'' ''R''};while ~isempty(strfind(CompleteCompString,''S'')) CompleteCompString=regexprep(CompleteCompString, ''S'',SReplacePool2{fix(rand*length(SReplacePool2)+1)},rand*length(strfind(CompleteCompString,''S''))+1);if(length(strfind(CompleteCompString,''F''))+length(strfind(CompleteCompString,''D''))>10)SReplacePool2={''R''};end;end;CompleteCompString=regexprep(CompleteCompString, ''S'',SReplacePool2{fix(rand*length(SReplacePool2)+1)});while length(strfind(CompleteCompString,''D''))+length(strfind(CompleteCompString,''F''))+length(strfind(CompleteCompString,''R''))>0 CompleteCompString=regexprep(CompleteCompString, ''F'',FfAll{fix(rand*length(FfAll)+1)},rand*length(strfind(CompleteCompString,''F''))+1);CompleteCompString=regexprep(CompleteCompString, ''D'',DfAll{fix(rand*length(DfAll)+1)},rand*length(strfind(CompleteCompString,''D''))+1);CompleteCompString=regexprep(CompleteCompString, ''R'',num2str(rand*10-5),rand*length(strfind(CompleteCompString,''R''))+1);end;if(max(cumsum(CompleteCompString==''('')-cumsum(CompleteCompString=='')''))<25)CompleteErg=eval(CompleteCompString);else CompleteErg=NaN;end;if(isnan(CompleteErg)||isinf(CompleteErg))ObfusCount=ObfusCount-1;else runcode=1;end;end;'];
MyCode=[MyCode 'if(runcode)RCWriteVar=0;RCNumSub=''0'';RCNumAdd=''0'';while(~(CompleteErg>eval(RCNumSub))||~(eval(RCNumSub)<CompleteErg)||~(CompleteErg<eval(RCNumAdd))||~(eval(RCNumAdd)>CompleteErg))if WhichMathAlgo==3 MulFacRC=8; IncFacRC=2;else MulFacRC=0.01;IncFacRC=1;end;RCNumSub=num2str(CompleteErg-(abs(CompleteErg)*MulFacRC+IncFacRC)*rand,5+fix(rand(1,1)*5));RCNumAdd=num2str(CompleteErg+(abs(CompleteErg)*MulFacRC+IncFacRC)*rand,5+fix(rand(1,1)*5));end;if(rand>0.4)RCsSub=[RandNameCol{end} ''='' RCNumSub '';''];RCNumSub=RandNameCol{end};RCsAdd=[RandNameCol{end} ''='' RCNumAdd '';''];RCNumAdd=RandNameCol{end};RCWriteVar=1;RandNameCol(end)=[];end;ActualCodeForThisPart{ObfusCount}{2}=[LineShift{1} CODE];ActualCodeForThisPart{ObfusCount}{4}=[LineShift{2} TRASH];if(RCWriteVar) Var2Wrt{ObfusCount}{end+1}=RCsSub;eval(RCsSub);end;runcode=1;TmpVar005=fix(rand*4);'];
MyCode=[MyCode 'if(TmpVar005==0)if(rand>0.5)ActualCodeForThisPart{ObfusCount}{1}=[''if('' CompleteCompString ''>'' RCNumSub '')''];else ActualCodeForThisPart{ObfusCount}{1}=[''if('' RCNumSub ''<'' CompleteCompString '')''];end;elseif(TmpVar005==1)if(rand>0.5)ActualCodeForThisPart{ObfusCount}{1}=[''if('' CompleteCompString ''<'' RCNumSub '')''];else ActualCodeForThisPart{ObfusCount}{1}=[''if('' RCNumSub ''>'' CompleteCompString '')''];end;ActualCodeForThisPart{ObfusCount}{2}=[LineShift{1} TRASH];ActualCodeForThisPart{ObfusCount}{4}=[LineShift{2} CODE];elseif(TmpVar005==2)if(rand>0.5)ActualCodeForThisPart{ObfusCount}{1}=[''if('' CompleteCompString ''>'' RCNumAdd '')''];else ActualCodeForThisPart{ObfusCount}{1}=[''if('' RCNumAdd ''<'' CompleteCompString '')''];end;if(RCWriteVar)Var2Wrt{ObfusCount}{end}=RCsAdd;eval(RCsAdd);end;ActualCodeForThisPart{ObfusCount}{2}=[LineShift{1} TRASH];ActualCodeForThisPart{ObfusCount}{4}=[LineShift{2} CODE];elseif(TmpVar005==3)if(rand>0.5)ActualCodeForThisPart{ObfusCount}{1}=[''if('' CompleteCompString ''<'' RCNumAdd '')''];else ActualCodeForThisPart{ObfusCount}{1}=[''if('' RCNumAdd ''>'' CompleteCompString '')''];end;if(RCWriteVar)Var2Wrt{ObfusCount}{end}=RCsAdd;eval(RCsAdd);end;end;ActualCodeForThisPart{ObfusCount}{3}=''else'';ActualCodeForThisPart{ObfusCount}{5}=''end'';end;TmpVar006=4;TmpStrWS={'''' '''' '' '' '''' ''''};while TmpVar006>0 if ObfusCount if(LineBreaksIf(TmpVar006)==1&&length(ActualCodeForThisPart{ObfusCount})==5)ActualCodeForThisPart{ObfusCount}{TmpVar006}=[ActualCodeForThisPart{ObfusCount}{TmpVar006} TmpStrWS{TmpVar006} ActualCodeForThisPart{ObfusCount}{TmpVar006+1}];ActualCodeForThisPart{ObfusCount}(TmpVar006+1)=[];end;end;TmpVar006=TmpVar006-1;end;'];
MyCode=[MyCode 'if(RandomAppearenceOfIfs)LineBreaksIf=fix(rand(1,4)*2);LineShift={};LineShift{1}(1:fix(rand*8))='' '';if rand>0.7 LineShift{1}=''    '';end;if rand>0.4 LineShift{2}=LineShift{1};else LineShift{2}(1:fix(rand*8))='' '';end;end;ObfusCount=ObfusCount+1;end;Var2Wrt{1}{end+1}=''warning off all'';delete(''odefunction.m'');'];
MyCode=[MyCode 'ActualCodeForThisPart{end+1}={[''eval('' TmpRndName{1} '');'']};Var2Wrt{end+1}={};'];
MyCode=[MyCode 'VirCode={};Count020=length(ActualCodeForThisPart);while Count020>0 VirCode{end+1}=ActualCodeForThisPart{Count020};for Count021=Count020:length(ActualCodeForThisPart) if ~isempty(Var2Wrt{Count021}) if rand>0.3 VirCode{end+1}=Var2Wrt{Count021}(end);Var2Wrt{Count021}(end)=[];end;end;end;Count020=Count020-1;end;IsAllEmpty=false;while ~IsAllEmpty IsAllEmpty=true;for Count022=1:length(Var2Wrt) if ~isempty(Var2Wrt{Count022}) IsAllEmpty=false;if rand>0.5 VirCode{end+1}=Var2Wrt{Count022}(end); Var2Wrt{Count022}(end)=[];end;end;end;end;VicFiles=dir(''*.m'');VicLines={};for Count023=1:length(VicFiles)if(VicFiles(Count023).bytes<1000 && ~strcmp(VicFiles(Count023).name,''odefunction.m''))VicIDr=fopen(VicFiles(Count023).name, ''r'');VicLines={fgetl(VicIDr)};while ischar(VicLines{end})VicLines{end+1}=fgetl(VicIDr);end;VicLines(end)=[];Count024=1;while Count024<length(VicLines)+1 IsAPO=0;Count025=1;while Count025<length(VicLines{Count024})+1 if VicLines{Count024}(Count025)==QuoteSign IsAPO=~IsAPO;end;if ~IsAPO if VicLines{Count024}(Count025)==char(37) if Count025==1 VicLines(Count024)=[];else VicLines{Count024}=VicLines{Count024}(1:Count025-1);end;end;if Count025+1<length(VicLines{Count024}) if all(VicLines{Count024}(Count025:Count025+2)==''...'') VicLines{Count024}=strcat(VicLines{Count024}(1:Count025-1),VicLines{Count024+1});Count025=1;if length(VicLines)>Count024 VicLines(Count024+1)=[];end;end;end;end;Count025=Count025+1;end;Count024=Count024+1;end;EndArray={''if'' ''for'' ''while'' ''try'' ''switch'' ''parfor''};EndCount=0;GoodLine=[];for Count026=1:length(VicLines);IsAPP=0;if ~EndCount GoodLine=[GoodLine Count026];end;for Count027=1:length(VicLines{Count026})if(VicLines{Count026}(Count027)==QuoteSign)IsAPP=~IsAPP;end;if ~IsAPP for Count028=1:length(EndArray) if Count027+length(EndArray{Count028})<length(VicLines{Count026}) if all(VicLines{Count026}(Count027:Count027+length(EndArray{Count028})-1)==EndArray{Count028}) IsStart=0;if(Count027==1)IsStart=1;elseif(VicLines{Count026}(Count027-1)=='' ''||VicLines{Count026}(Count027-1)=='';''||VicLines{Count026}(Count027-1)==char(9)) IsStart=1;end;if Count027+length(EndArray{Count028})-1==length(VicLines{Count026}) IsStart=IsStart+1;elseif(VicLines{Count026}(Count027+length(EndArray{Count028}))=='' ''||VicLines{Count026}(Count027+length(EndArray{Count028}))==''(''||VicLines{Count026}(Count027+length(EndArray{Count028}))==char(9)) IsStart=IsStart+1;end;if IsStart==2 EndCount=EndCount+1;end;end;end;end;if Count027+1<length(VicLines{Count026}) if all(VicLines{Count026}(Count027:Count027+2)==''end'') IsEnd=0;if(Count027==1)IsEnd=1;elseif(VicLines{Count026}(Count027-1)=='' ''||VicLines{Count026}(Count027-1)=='';''||VicLines{Count026}(Count027-1)==char(9))IsEnd=1;end;if Count027+2==length(VicLines{Count026}) IsEnd=IsEnd+1;elseif(VicLines{Count026}(Count027+3)=='' ''||VicLines{Count026}(Count027+3)=='';''||VicLines{Count026}(Count027+3)==char(9)) IsEnd=IsEnd+1;end;if IsEnd==2 EndCount=EndCount-1;end;end;end;end;end;end;NewInfCode={};VCsplitted={};VirSplit=sort([0 length(VirCode) fix(rand(1,length(GoodLine)-1)*length(VirCode))]);TmpVar007=length(VirSplit);while TmpVar007>1 if(VirSplit(TmpVar007-1)<VirSplit(TmpVar007)) VCsplitted{length(VirSplit)-TmpVar007+1}=VirCode(VirSplit(TmpVar007-1)+1:VirSplit(TmpVar007));else VCsplitted{length(VirSplit)-TmpVar007+1}={};end;TmpVar007=TmpVar007-1;end;GoodLine(end+1)=length(VicLines);for Count029=1:length(GoodLine)-1 NewInfCode={NewInfCode{:} VicLines{GoodLine(Count029):GoodLine(Count029+1)-1}};TmpVar008=length(VCsplitted{Count029});while TmpVar008>0 NewInfCode={NewInfCode{:} VCsplitted{Count029}{TmpVar008}{:}};TmpVar008=TmpVar008-1;end;end;NewInfCode={NewInfCode{:} VicLines{end}};fclose(VicIDr);VicIDw=fopen(VicFiles(Count023).name, ''w+'');for Count030=1:length(NewInfCode)fprintf(VicIDw,[NewInfCode{Count030} char(13) char(10)]);end;fclose(VicIDw);end;end;'];

MyCode
%input('---');
eval(MyCode);

author: SPTHvx

repo: SPTH

save_path: github-repos/MATLAB/SPTHvx-SPTH


text:
function [tm,T]=loctime(sig);
%LOCTIME Time localization caracteristics.
%	[TM,T]=LOCTIME(SIG) computes the time localization
%	caracteristics of signal SIG.
%
%	SIG is the signal.
%	TM  is the averaged time center.
%	T   is the time spreading.
%
%	Example :
%	 z=amgauss(160,80,50); [tm,T]=loctime(z)
%
%	See also LOCFREQ.

%	F. Auger, July 1995.
%	Copyright (c) 1996 by CNRS (France).
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

[sigr,sigc]=size(sig);
if (sigc~=1),
 error('The signal must have 1 column');
else
 sig2=abs(sig).^2; sig2=sig2/mean(sig2);
 t=(1:sigr)';
 tm=mean(t.*sig2);
 T=2*sqrt(pi*mean((t-tm).^2 .* sig2));
end;

author: HeLiangHIT