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

    methods
        function obj = Tx(varargin)
            obj@ETS3('Tx', varargin{:});

        end

        % function s = char(ets)
        %     s = ets.what;
        % end
    end
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function normal_ms_cdf_test ( )

%*****************************************************************************80
%
%% NORMAL_MS_CDF_TEST tests NORMAL_MS_CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 March 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NORMAL_MS_CDF_TEST\n' );
  fprintf ( 1, '  NORMAL_MS_CDF evaluates the CDF\n' );
  fprintf ( 1, '  for the Normal MS distribution.\n' );

  mu = 100.0;
  sigma = 15.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter MU =    %14f\n', mu );
  fprintf ( 1, '  PDF parameter SIGMA = %14f\n', sigma );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '       X            CDF      CDF_INV\n' );
  fprintf ( 1, '\n' );

  for i = -20 : 20
    x = mu + sigma * i / 10.0;
    cdf = normal_ms_cdf ( x, mu, sigma );
    x2 = normal_ms_cdf_inv ( cdf, mu, sigma );
    fprintf ( 1, '  %14g  %14g  %14g\n', x, cdf, x2 );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Fig. 8.5   Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
%

clear all;
close all;

th=0:.1:pi;
Im=sin(th);
Re=cos(th);

figure(1)
plot(Re,Im,'-',.5*Re,.5*Im,'--'),grid
axis([-1 1 -.5 1.5])
axis('square')
hold on
plot(1,0,'x')
plot(.6,0,'x')
plot(.3,0,'x')
plot(0,0,'x')
plot(-.5,0,'x')
plot(-1,0,'x')
j=sqrt(-1);
p=[exp(j*pi/4) exp(-j*pi/4)];
plot(p,'x')
p2=[exp(j*pi/2) exp(-j*pi/2)];
plot(p2,'x')
p3=1.2*[exp(j*pi/2) exp(-j*pi/2)];
plot(p3,'x')
p4=.5*[exp(j*pi/4) exp(-j*pi/4)];
plot(p4,'x')
p5=.5*[exp(j*pi/2) exp(-j*pi/2)];
plot(p5,'x')
p6=.5*[exp(j*3*pi/4) exp(-j*3*pi/4)];
plot(p6,'x')
title('Fig. 8.5  Time sequences associated with poles in the z-plane')
hold off

figure(2) % This figure creates four of the time responses
num=[1  0];
den=[1 -1];
N=11;
n=0:N-1;
y=dimpulse(num,den,N);
Ny=size(y);
subplot(2,2,1)
axis([1 10 -1.5 1.5])
zero=zeros(Ny);
plot(n,y,'*',n,y,'--',n,zero,'-')
title('z = +1')

den=[1 -.6];
subplot(2,2,2)
y=dimpulse(num,den,N);
plot(n,y,'*',n,y,'--',n,zero,'-')
title('z = +.6')

den=[1 -.3];
subplot(2,2,3)
y=dimpulse(num,den,N);
plot(n,y,'*',n,y,'--',n,zero,'-')
title('z = +.3')

den=[1 0];
subplot(2,2,4)
y=dimpulse(num,den,N);
plot(n,y,'*',n,y,'--',n,zero,'-')
title('z = 0')

figure(3)   % This figure creates four more of the time responses
num=[1  0];
j=sqrt(-1);
p=[exp(j*pi/4) exp(-j*pi/4)];
den=poly(p);
N=11;
n=0:N-1;
y=.67*dimpulse(num,den,N);
Ny=size(y);
subplot(2,2,1)
axis([1 10 -1.5 1.5])
zero=zeros(Ny);
plot(n,y,'*',n,y,'--',n,zero,'-')
title('z = exp(+-j*45 deg)')

p=[exp(j*pi/2) exp(-j*pi/2)];
den=poly(p);
subplot(2,2,2)
y=dimpulse(num,den,N);
plot(n,y,'*',n,y,'--',n,zero,'-')
title('z = exp(+-j*90 deg)')

den=[1  1];
subplot(2,2,3)
y=dimpulse(num,den,N);
plot(n,y,'*',n,y,'--',n,zero,'-')
title('z = -1')

p=1.2*[exp(j*pi/2) exp(-j*pi/2)];
den=poly(p);
subplot(2,2,4)
y=.5*dimpulse(num,den,N);
plot(n,y,'*',n,y,'--',n,zero,'-')
title('z = 1.2*exp(+-j*90 deg)')

figure(4)  % This figure creates four more of the time responses
num=[1  0];
j=sqrt(-1);
p=.5*[exp(j*pi/4) exp(-j*pi/4)];
den=poly(p);
N=11;
n=0:N-1;
y=dimpulse(num,den,N);
Ny=size(y);
subplot(2,2,1)
axis([1 10 -1.5 1.5])
zero=zeros(Ny);
plot(n,y,'*',n,y,'--',n,zero,'-')
title('z = .5*exp(+-j*45 deg)')

p=.5*[exp(j*pi/2) exp(-j*pi/2)];
den=poly(p);
subplot(2,2,2)
y=dimpulse(num,den,N);
plot(n,y,'*',n,y,'--',n,zero,'-')
title('z = .5*exp(+-j*90 deg)')

p=.5*[exp(j*3*pi/4) exp(-j*3*pi/4)];
den=poly(p);
subplot(2,2,3)
y=dimpulse(num,den,N);
plot(n,y,'*',n,y,'--',n,zero,'-')
title('z = .5*exp(+-j*135 deg)')

den=[1 .5];
subplot(2,2,4)
y=dimpulse(num,den,N);
plot(n,y,'*',n,y,'--',n,zero,'-')
title('z = -.5')

author: Sable

repo: mcbench-benchmarks


text:
%helpme_ch4    helpme  for Chapter 4 specialized functions
%  IFISS scriptfile: HCE; 27 April 2005.

fprintf(' \n');
fprintf(' The files in the directory solvers/ch4_code can be used to reproduce \n');
fprintf(' results from Chapter 4 of the Elman-Silvester-Wathen monograph.\n');
fprintf(' Type "help ch4_code" for details.\n');
fprintf(' \n');


text:
function varargout=lpcrand(varargin)
%
% For calling details please see v_lpcrand.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_lpcrand.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: lpcrand.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_lpcrand(varargin{:});
else
    v_lpcrand(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
clear; clc;
load mariodata;

global H;
global A;
global mario;
global dxnkey;
global dxnpressed;
global shiftpressed;
global spacepressed;
global spacetime;
dxnpressed=0;
shiftpressed=0;
spacepressed=0;
dxnkey='';

timeout	= 30;	% duration of sim (sec)
% rfrsh	= 1/60;	% NTSC refresh rate (sec)
rfrsh	= 1/10;	% refresh rate that works for now (want better rate)
scale	= 1;	%

%SETUP FIGURE:
figure(1); clf;
hold on;
gnd=patch([0 0 10000 10000 0],[-20 0 0 -20 -20],7);
set(gca,'YTick',[],'XTick',[],'Color',[0.3 0.6 0.94]);
plot3([0 1000],[0 0],[-1 -1],'k','LineWidth',2);
axis([0 400 -20 200]);
set(gcf,'keypressfcn',@keypress,'keyreleasefcn',@keyrelease,'WindowButtonUpFcn','');
axis equal;

%INITIALIZE SPRITE PATCHES:
maxdims = [32 17]; %sprite dimensions
colormap(clrs);
for i=1:maxdims(1)	%rows
	for j=1:maxdims(2)	%cols
		H(i,j)=patch([-1 -1 0 0 -1]+j,[-1 0 0 -1 -1]+i,0,'FaceAlpha',0);
	end;
end;
shading flat;
ttl = title(sprintf('TIME\n%03.0f',max(timeout,0)),'FontWeight','bold','Color',[1 1 1],'FontSize',14,'FontName','Monotxt','HorizontalAlignment','right','Position', [390,150,1]);

%INITIAL STATE:
%states:	1:small, 2:big, 3:firepower
%actions:	1:stand, 2:run, 3:jump, 4:duck, 5:skid, 6:climb, 7:switm, 8:shoot, 9: die
%dxns:		1:right, 2:left
state	= 3;		%current character state
dxn		= 1;		%current character direction
action	= 1;		%current character action
iter	= 1;		%frame within current action
loc		= [0 0];	%initial position
spd		= [0 0];	%current x and y speed
spdlim	= [8 18];	%magnitude of max speed in x and y
acc		= [0 0];	%current x and y acceleration
updatesprite(state,dxn,action,iter,[0 0])

%REMINDER:
%states:	1:small	2:big  3:firepower
%dxns:		1:right 2:left
%actions:	1:stand 2:walk 3:jump 4:duck 5:skid 6:climb 7:swim 8:shoot 9:die

%MAIN LOOP:
tic;
tnext=0;
t=toc;
while toc<timeout,
	%DIRECTION COMMANDED:
	if dxnpressed,				%accelerate
		if dxnkey(1)=='d'	%user pressed down (duck)
			action=4;
			if spd(1)==0	%stay at rest
				acc(1)=0;
			else
				acc(1) = -sign(spd(1))*ceil(abs(spd(1)/2)); %decellerate
			end;

		elseif dxnkey(1)=='r'		%accelerate right
			acc(1) = 2*(spd(1)<spdlim(1));
			if spd(1)<0			%skidding turn left to right
				action	= 5;
				dxn		= 2;
			else				%running right
				action	= 2;
				dxn		= 1;
			end;

		elseif dxnkey(1)=='l'	%accelerate left
			acc(1) = -2*(spd(1)>-spdlim(1));
			if spd(1)>0			%skidding turn right to left
				action	= 5;
				dxn		= 1;
			else				%running left
				action	= 2;
				dxn		= 2;
			end;
		end;
	%NO ARROWS PRESSED: slow to rest
	else
		if spd(1)==0	%stay at rest
			acc(1)=0;
			action=1;
		else
			acc(1) = -sign(spd(1))*ceil(abs(spd(1)/2)); %decellerate
		end;
	end;

	%SHOOTING:
	if state==3 && shiftpressed,
		action=8;
	end;

	%Jumping:
	if spacepressed && loc(2)==0,
		spacepressed = 0;
		acc(2)=-3;
		spd(2)=18;
	end;

	%UPDATE POSITION:
	spd = round(spd+acc);
	loc = loc+spd;

	%end jump?
	if loc(2)<=0 && spd(2)<0	%end jump
		loc(2)=0;
		spd(2)=0;
		acc(2)=0;
		if dxnkey(1)=='l' || dxnkey(1)=='r' %end running
			action=2;
		elseif dxnkey(1)=='d'				%end ducking
			action=4;
		end;
	elseif loc(2)>0
		if action==4 || action==8 %jumping stance unless ducking or shooting
		else
			action=3;	%airborne stance
		end;
	end;

	%INCREMENT FRAME OF CURRENT ACTION:
	if action==2,
		if iter==1, iter=3;
		else iter=iter-1;
		end;
	elseif sum(action==[1 3 4 5 7 9])
		iter=1;
	end;

	%UPDATE FIGURE:
	updatesprite(state,dxn,action,iter,loc)
	axis([0 400 -20 200])
	set(ttl,'String',sprintf('TIME\n%03.0f',max(timeout-toc,0)));
	drawnow;

	%CONTROL REFRESH RATE:
% 	tnext	= tnext+rfrsh;
% 	t=toc;
% 	txtra	= tnext-t;
% 	if txtra>0,
% 		pause(txtra);
% 		disp('paused');
% 	else
% 		disp([t,tnext,txtra])
% 	end;
end;

%DEATH
ylims = get(gca,'YLim');
spd=[0 12];
acc=[0 -3];
updatesprite(1,1,9,1,loc);
pause(1);
while loc(2)>ylims(1)-30
	spd(2)=max(spd(2)+acc(2),-spdlim(2)); %terminal velocity
	loc=loc+spd;
	updatesprite(1,1,9,1,loc);
	drawnow;
end;

clf;
plot(0,0);
set(gca,'Color','k','XTick',[],'YTick',[],'ButtonDownFcn','test1');
text(140,100,'click to continue','Color','w','FontSize',18);
axis([0 400 -20 200],'equal');


text:
% CELLSRC2DB Constitute a database compatible with the ScatNet
%   convention with a cell of one object per image and a src
%
% Usage
%   db = cellsrc2db(cell, src);
function db = cellsrc2db(cell, src)
	db.src = src;
	db.features = cell2mat(cell);
	for i = 1:numel(cell)
		db.indices{i} = i;
	end
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function test_gt(prediction_folder, view_label_folder, write_result)
% prediction_folder contains <classname>_pred_view.txt
% view_label_folder contains <classname>.txt
% write_result - if set to 1, write result to acc_mederr_results.txt in prediction_folder

addpath(fullfile(mfilename('fullpath'), '../../'));
global_variables;

N = numel(g_cls_names);

acc = zeros(1,N);
mederr = zeros(1,N);

for k = 1:N
    cls = g_cls_names{k};
    try
        [acc(k), mederr(k)] = compute_vp_acc_mederror(fullfile(prediction_folder,sprintf('%s_pred_view.txt', cls)), ...
            fullfile(view_label_folder, sprintf('%s.txt', cls)));
    catch
    end
end
display(acc);
display(mederr);
fprintf('Mean acc = %f\n', mean(acc));
fprintf('Mean mederr = %f\n (in degree)\n', mean(mederr)/pi*180);

if nargin < 3
    write_result = 1;
end
if write_result
    fid = fopen(fullfile(prediction_folder, 'acc_mederr_results.txt'), 'w');
    for k = 1:N
        fprintf(fid, sprintf('%s ', g_cls_names{k}));
    end
    fprintf(fid, '\n');
    for k = 1:N
        fprintf(fid, sprintf('%f ', acc(k)));
    end
    fprintf(fid, '\n');
    for k = 1:N
        fprintf(fid, sprintf('%f ', mederr(k)/pi*180));
    end
    fprintf(fid, '\n');
    fprintf(fid, 'Mean acc = %f\n', mean(acc));
    fprintf(fid, 'Mean mederr = %f\n (in degree)\n', mean(mederr)/pi*180);
    fclose(fid);
end

author: ShapeNet

repo: RenderForCNN


text:
% easy11    Creates a stereographic plot of GPS satellite  orbits
%                from an epehemerides or almanac file. The plot is as
%                seen from the position (phi, lambda). All (parts of)
%                orbits with elevation angle lower than a given
%                cut-off angle (mask) are omitted.
%                As input we use the ephemerides in our basic RINEX
%               navigation file. It only contains nine ephemerides. The
%               result shows two windows with a considerable amount
%               of visible satellites. This is like the actual situation in the
%               early 1990es. Then you could have a window  of about
%               one hour at say 5--6 am and this repeated 12 hours later.
%
%              An additional plot is created showing number
%              of visible satellites and when they are visible.
%              Finally a stem plot depicts the number of
%             satelliets visible for how long
%
%             Alternatively, if you want to include all GPS satellites,
%             we need an almanac file which can be downloaded from
%             http://www.navcen.uscg.gov/?pageName=gpsAlmanacs

%Kai Borre 26-08-2008
%Copyright (c) by Kai Borre
%$Revision: 1.1 $  $Date: 2009/04/05  $
% Total revision January 10, 2016

% RINEX version 3.03

set(0,'DefaultTextFontName','Times');
set(0,'DefaultAxesFontName','Times');
set(0,'DefaultTextFontSize',12);

% If you want to use a RINEX ephemeris file uncomment
% line 37 and comment line 41 out.

% the following five assignments are mandatory
%%rinexe('log_24h.15n','eph.dat');
% in  case you want to include all GPS satellites in the constellation,
% we need an almanac. This is read and converted to the usual
% eph format by the following call
rinexa('current.alm','eph.dat');

ephemerides = 'eph.dat';
mask = 10;
phi = [57 12 43];
lambda = [50 10 39];

%reading ephemerides
fide = fopen(ephemerides,'r');
Eph = fread(fide,inf,'double');
m = length(Eph);
eph = reshape(Eph,21,m/21);

% transformation of given location (phi,lambda,h) to (X,Y,Z)
Phi = dms2rad(phi(1),phi(2),phi(3));
Phi = Phi*180/pi;
Lambda = dms2rad(lambda(1),lambda(2),lambda(3));
Lambda = Lambda*180/pi;
[M(1,1),M(2,1),M(3,1)] = frgeod(6378137,298.257222101,Phi,Lambda,0);

% Computation of (azimuth, elevation) for each satellite
% for each 15 min.s. We use only one ephemeris for each PRN.
% Anyway, the plot is only for visual use
[prns, ind] = unique(eph(1,:));
satrows = length(ind);
 az = ones(satrows,96)*inf; % satrows is number of PRN and 96 quaters of an hour in 24 hours
el = ones(satrows,96)*inf;

start_time = max(eph(21,:));
for sat = ind'
    j = 0;
    for time = start_time:900:start_time+86400
        S = satpos(time,eph(:,sat)); %sat
        j = j+1;
        [azimuth,elevation,distance] = topocent(M,S-M);
        az(sat,j) = azimuth; % sat runs over PRN, j runs over sow
        el(sat,j) = elevation;
    end
end
% We assume that there are satrows PRNs
XX = zeros(satrows,40)*inf; % a matrix of NaNs to store plot data
YY = XX;

figure(1);
% polarhg draws coordinate lines of a polar plot. We add
% circles with radii 30 and 60 degrees
polarhg([30 60])
hold on
for k = ind' % 1:32
    % if az(k,1) == 0, break, end  % continue
    AZ = az(k,:);
    EL = el(k,:);
    % remove data below the cut-off angle
    AZ(find(EL <= mask)) = nan;
    EL(find(EL <= mask)) = nan;
    % conversion from polar to rectangular coordinates
    xx = (90-EL).*cos(AZ*pi/180);
    yy = (90-EL).*sin(AZ*pi/180);
    XX(k,1:length(xx)) = xx;
    YY(k,1:length(yy)) = yy;
end % k
%the first coord. moves text vertically (increasing values up),
% the second coord. moves text horizontally (increasing values right)
text(135,-95,{['Skyplot for the position (\phi, \lambda) = (' ...
    num2str(round(Phi)) '\circ, '  num2str(round(Lambda)) '\circ)']})
text(115,-45,{['Elevation mask  ' num2str(mask) '\circ' ]}) %120
text(-120,-120,['All PRNs except  ' num2str(setdiff(1:32,prns)) ])
plot(XX',YY','linewidth',2)
hold off
set(gca,'Fontsize',9);

print('easy111','-dpdf')

% preparation for visibility plot  %%%

% we choose a resolution of 5 min.s,
% ie. 24 hours times 12 = 288 which becomes the range of j
satsum = zeros(1,288);
visible = zeros(2*(size(prns,2)+1),288);

for sat = ind'
    %Az = [];
    %El = [];
    i = eph(1,sat);
    for j = 1:288
        time = 300*(j-1); % in s, 0<time<86400
        S = satpos(time,eph(:,sat));
        [az,el,d] = topocent(M,S-M);
        if el > mask
 %           Az = [Az az];
  %          El = [El el];
            satsum(j) = satsum(j)+1;
            visible(2*i,j) = 1;
        end
    end
end

figure(2);
set(gca,'Fontsize',16);
area(satsum)
set(gca,'XTick',1:71:288)
set(gca,'XTickLabel',{'0','6','12','18','24'})
xlabel('GPS Time [hours]')
ylabel('# of Visible Satellites')
title(['Elevation Mask ' num2str(mask) '\circ'])
colormap summer

print -dpdf easy112

figure(3);
set(gca,'Fontsize',16);
imagesc(flipud(visible));
colormap(autumn(5))
set(gca,'XTick',1:71:288)
set(gca,'XTickLabel',{'0','6','12','18','24'})
set(gca,'YTick',-3:16:(2*(size(prns,2)+1)))
set(gca,'YTickLabel',{'2','8','16','24','32'});
xlabel('GPS Time [hours]')
ylabel('PRNs')
title('Yellow Lines Indicate Visible Satellites')
colormap winter

print -dpdf easy113

figure(4);
set(gca,'Fontsize',16);
an = sum(visible,1);
% we assume the total # of PRNs always is between 1 and 20
total = zeros(1,20);
for i = 1:288
    j = an(i);
    if j == 0, j = j+1; total(j) = 0.01; continue, end
    total(j) = total(j)+1;
end
% the adder total adds per 5 minutes, division by 12 changes to per hour
bar(1:20,total/12);  % 10-> 32
title('Number of Visible Satellites')
ylabel('Hours')

print -dpdf easy114

%%%%%%%%%%%%%%%%%%%%% end easy11.m %%%%%%%%%%%%%%

author: yandld


text:
function [KFSmatrices] = makeMatricesSLMfullSV(dataValues,p)
%% Purpose
%  This function calculates the required matrices for the Kalman Filter for the SLM full SV model.
% The model is Y

%% Initialize
Ppsi=dataValues.Ppsi;
[ns,n]=size(Ppsi);

%% Z: coefficients of measurment equation
Ztemp=zeros(ns+n,n*(p+1));
Ztemp(ns+1:end,1:2*n)=[eye(n) eye(n)];
Ztemp(1:ns,1:n)= Ppsi;
KFSmatrices.Zmatrix=Ztemp;

%% R: for variance covariance in transition equation
KFSmatrices.R=[eye(2*n);zeros(n*(p-1),2*n)];

%% L: for variance covariance in measurement equation
KFSmatrices.L=[eye(ns);zeros(n,ns)];

%% Tmatrix: coefficients in transition equation
TmatrixTemp=zeros(n*(p+1),n*(p+1));
TmatrixTemp(1:n,1:n) = eye(n);
%TmatrixTemp(n+1:2*n,n+1:end) = Bhat';
TmatrixTemp(2*n+1:end,n+1:end-n) = eye(n*(p-1));

KFSmatrices.Tmatrix=TmatrixTemp;

end

author: european-central-bank

repo: BEAR-toolbox


text:
%                               getExtraOpts.m
%
% This function creates and outputs a struct that consists of all the
% options in opts but not in otherOpts.
% It is used as a helper function in manageOptions.m and
% manageOptionsForBenchmark.m
%
% Inputs:
%         opts(struct)       :  consists of options.
%         otherOpts(struct)  :  consists of options.
% Outputs:
%         extras(struct)     :  consists of extral options appearing in
%                               opts but not in otherOpts.
%
%
% PhasePack by Rohan Chandra, Ziyuan Zhong, Justin Hontz, Val McCulloch,
% Christoph Studer, & Tom Goldstein
% Copyright (c) University of Maryland, 2017

function extras = getExtraOpts(opts, otherOpts)
    extras = struct;
    optNames = fieldnames(opts);

    for i = 1 : length(optNames)
        optName = optNames{i};
        if ~isfield(otherOpts, optName)
            extras.(optName) = opts.(optName);
        end
    end
end

author: tomgoldstein

repo: phasepack-matlab


text:
classdef DesignVarMonitor_LevelSet < DesignVarMonitor_Abstract

    properties (Access = protected, Abstract)
        unfittedType
        meshIncludeBoxContour
    end

    properties (Access = protected)
        designVarName = 'Level Set - \phi';
        meshUnfitted
    end

    methods (Access = public)

        function obj = DesignVarMonitor_LevelSet(cParams)
            obj@DesignVarMonitor_Abstract(cParams);
        end

        function plot(obj)
            obj.refreshFigure();
        end

    end

    methods (Access = protected)

        function initPlotting(obj)
             set(obj.axes,'CLim',[0, 1],'XTick',[],'YTick',[]);
        end

    end

    methods (Access = protected, Static)

        function color = getColor()
            color = [1 0 0];
        end

    end

    methods (Access = private)

        function refreshFigure(obj)
            figure(obj.figHandle.Number)
            cla reset;
            hold on
            uMesh = obj.designVar.getUnfittedMesh;
            uMesh.plot();
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%[2018]-"Butterfly optimization algorithm: a novel approach for global
%optimization"

% (9/12/2020)

function BOA = jButterflyOptimizationAlgorithm(feat,label,opts)
% Parameters
lb    = 0;
ub    = 1;
thres = 0.5;
c     = 0.01;   % modular modality
p     = 0.8;    % switch probability

if isfield(opts,'T'), max_Iter = opts.T; end
if isfield(opts,'N'), N = opts.N; end
if isfield(opts,'c'), c = opts.c; end
if isfield(opts,'p'), p = opts.p; end
if isfield(opts,'thres'), thres = opts.thres; end

% Objective function
fun = @jFitnessFunction;
% Number of dimensions
dim = size(feat,2);
% Initial
X   = zeros(N,dim);
for i = 1:N
	for d = 1:dim
    X(i,d) = lb + (ub - lb) * rand();
	end
end
% Pre
Xnew = zeros(N,dim);
fitG = inf;
fit  = zeros(1,N);

curve = inf;
t = 1;
% Iterations
while t <= max_Iter
  % Fitness
  for i = 1:N
    fit(i) = fun(feat,label,(X(i,:) > thres),opts);
    % Global update
    if fit(i) < fitG
      fitG = fit(i);
      Xgb  = X(i,:);
    end
  end
  % Power component, increase from 0.1 to 0.3
  a = 0.1 + 0.2 * (t / max_Iter);
  for i = 1:N
    % Compute fragrance (1)
    f = c * (fit(i) ^ a);
    % Random number in [0,1]
    r = rand();
    if r < p
      r1 = rand();
      for d = 1:dim
        % Move toward best butterfly (2)
        Xnew(i,d) = X(i,d) + ((r1 ^ 2) * Xgb(d) - X(i,d)) * f;
      end
    else
      % Random select two butterfly
      R  = randperm(N);
      J  = R(1);
      K  = R(2);
      r2 = rand();
      for d = 1:dim
        % Move randomly (3)
        Xnew(i,d) = X(i,d) + ((r2 ^ 2) * X(J,d) - X(K,d)) * f;
      end
    end
    % Boundary
    XB = Xnew(i,:); XB(XB > ub) = ub; XB(XB < lb) = lb;
    Xnew(i,:) = XB;
  end
  % Replace
  X = Xnew;
  % Save
  curve(t) = fitG;
  fprintf('\nIteration %d Best (BOA)= %f',t,curve(t))
  t = t + 1;
end
% Select features
Pos   = 1:dim;
Sf    = Pos((Xgb > thres) == 1);
sFeat = feat(:,Sf);
% Store results
BOA.sf = Sf;
BOA.ff = sFeat;
BOA.nf = length(Sf);
BOA.c  = curve;
BOA.f  = feat;
BOA.l  = label;
end

author: JingweiToo

repo: Wrapper-Feature-Selection-Toolbox


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
% Problem 4
% Inverse z-Transform computation of X(z)=(2z+3)/(z^2+5z+6)

% a)in rational form
syms n z
X=(2*z+3)/(z^2+5*z+6);
iztrans(X,n)

% b) in partial fraction form
num=[ 2 3]
den= [ 1 5 6];
[R,P,K]=residue(num,den)

X=3/(z+3)-1/(z+2)
iztrans(X,n)


text:
function [simulatedata,ht,pseudorc]=matrix_garch_simulate(t,k,parameters,p,o,q,m)
% Simulation of symmetric and asymmetric MATRIX multivariate GARCH models
%
% USAGE:
%   [SIMULATEDATA, HT, PSEUDORC] = matrix_garch_simulate(T, K, PARAMETERS, P, O, Q, M)
%
% INPUTS:
%   T            - Length of the time series to be simulated
%   K            - Cross-sectional dimension
%   PARAMETERS   - A 3-D matrix of K by K matrices where
%                    CC' = PARAMETERS(:,:,1), AA'(j)=PARAMETERS(:,:,1+j),
%                    GG'(j) = PARAMETERS(:,:,1+P+j), BB'(j)=PARAMETERS(:,:,1+P+Q+j)
%                    -OR- a K(K+1)/2*(1+P+O+Q) by 1 vector of parameters of
%                    the form returned my calling matrix_garch
%   P            - Positive, scalar integer representing the number of lags of the innovation process
%   O            - Non-negative scalar integer representing the number of asymmetric lags to include
%   Q            - Non-negative scalar integer representing the number of lags of conditional covariance
%   M            - [OPTIONAL] Number of ``intradaily'' returns to simulate to pseudo-Realized
%                    Covariance. If omitted, set to 72.
%
% OUTPUTS:
%   SIMULATEDATA - A time series with constant conditional correlation covariance
%   HT           - A [k k t] matrix of simulated conditional covariances
%   PSEUDORC     - A [k k t] matrix of pseudo-Realized Covariances
%
% COMMENTS:
%    The conditional variance, H(t), of a MATRIX GARCH is modeled as follows:
%
%      H(t) = CC' + AA'(1).*r_{t-1}'*r_{t-1} + ... + AA'(P).*r_{t-P}'*r_{t-P}
%                 + GG(1)'.*n_{t-1}'*n_{t-1} + ... + GG(O)'.*n_{t-P}'*n_{t-P}
%                  + BB(1)'.*H(t-1) +...+ BB(Q)'.*H(t-q)
%
%    where n_{t} = r_{t} .* (r_{t}<0).  If using realized measures, the
%    RM_{t-1} replaces r_{t-1}'*r_{t-1}, and the asymmetric version
%    replaces n_{t-1}'*n_{t-1}
%
%   Pseudo Realized Covariances are simulated by generating m-intra daily returns from a N(0,1/m)
%   and computing the Realized Covariance of these. These were used in Patton and Sheppard (2009)
%   when evaluating variance and covariance specifications in a Monte Carlo.  If M=1, then PSEUDORC
%   is just the outer product of the SIMULATEDATA.
%
%   NOTE: This program generates 2000 more than required to minimize any start-up bias
%
% See also TARCH_SIMULATE, CCC_GARCH_SIMULATE, MATRIX_GARCH

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch nargin
    case 6
        m=[];
    case 7
        % nothing
    otherwise
        error('MFEToolbox:Input','6 or 7 inputs required.');
end

% t
if ~isscalar(t) || t<0 || floor(t)~=t
    error('T must be a positive integer.');
end
% k
if ~isscalar(k) || k<2 || floor(k)~=k
    error('K must be a positive integer greater than or equal to 2.');
end
% p
if ~isscalar(p) || p<1 || floor(p)~=p
    error('P must be a positive integer if scalar.');
end

% o

if ~isscalar(o) || o<0 || floor(o)~=o
    error('O must be a non-neagative integer if scalar.');
end
% q

if ~isscalar(q) || q<0 || floor(q)~=q
    error('Q must be a non-neagative integer if scalar.');
end
% m
if isempty(m)
    m = 72;
end
if ~isscalar(m) || floor(m)~=m || m<1
    error('M must be a positive integer.')
end
% parameters
k2 = k*(k+1)/2;
if ismatrix(parameters)
    parameterCount = k*(k+1)/2 *(1+ p + o + q);
    if size(parameters,2)>size(parameters,1)
        parameters = parameters';
    end
    if length(parameters)~=parameterCount
        error('PARAMETERS must be K(K+1)/28(1+P+O+Q) when using a vector.')
    end
    parameterMatrices= zeros(k,k,1+p+o+q);
    index = 0;
    for i=1:(1+p+o+1)
        temp = vec2chol(parameters(index+1:index+k2));
        parameterMatrices(:,:,i) = (temp*temp'+temp*temp')/2;
        index=index+k2;
    end
    parameters = parameterMatrices;
elseif ndims(parameters)==3
    if any(size(parameters)~=[k k 1+p+o+q])
        error('PARAMETERS must be K by K by (1+P+O+Q) when using a 3-D matrix.')
    end
    parameterMatrices = parameters;
else
    error('The size of PARAMETERS is not compatible with this function.')
end

% Check stationarity
sumParameters = zeros(k);
for i=2:size(parameterMatrices,3)
    if i<=(p+1) || i>(1+p+o)
        w=1;
    else
        w=0.5;
    end
    sumParameters = sumParameters + w*parameterMatrices(:,:,i);
end
sumParameters = (sumParameters+sumParameters')/2;
if max(diag(sumParameters))>1
    warning('MFE::Nonstationary','The parameters do not correspond to a stationary solution.  Check for overflow.')
    uncond = parameters(:,:,1)/.005;
else
    uncond = parameters(:,:,1)./(ones(k)-sumParameters);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Set the burnin amount, 2000 is probably reasonable
burnin=2000;
%Up t by the burnin amount
t=t+burnin;

%Draw some normal random numbers
intraRandomNums=randn(m*t,k)*sqrt(1/m);
randomNums = cumsum(intraRandomNums);
randomNums = diff([zeros(1,k);randomNums(m:m:m*t,:)]);

%Perform the recursion
backCast = uncond;
backCastAsym = uncond;
T = t;
Ht = zeros(k,k,T);
simulatedData = zeros(T,k);
for t=1:T;
    Ht(:,:,t)=parameterMatrices(:,:,1);
    for j=1:p
        if (t-j)<1
            Ht(:,:,t)=Ht(:,:,t)+parameterMatrices(:,:,j+1).*backCast;
        else
            r = simulatedData(t-j,:);
            Ht(:,:,t)=Ht(:,:,t)+parameterMatrices(:,:,j+1).*(r'*r);
        end
    end
    for j=1:o
        if (t-j)<1
            Ht(:,:,t)=Ht(:,:,t)+parameterMatrices(:,:,p+j+1).*backCastAsym;
        else
            n = simulatedData(t-j,:).*(simulatedData(t-j,:)<0);
            Ht(:,:,t)=Ht(:,:,t)+parameterMatrices(:,:,p+j+1).*(n'*n);
        end
    end
    for j=1:q
        if (t-j)<1
            Ht(:,:,t)=Ht(:,:,t)+parameterMatrices(:,:,p+o+j+1).*backCast;
        else
            Ht(:,:,t)=Ht(:,:,t)+parameterMatrices(:,:,p+o+j+1).*Ht(:,:,t-j);
        end
    end
    simulatedData(t,:) = randomNums(t,:)*Ht(:,:,t)^(0.5);
end

%Initialize the covariance
pseudorc=zeros(k,k,t);
for t=1:T
    r = intraRandomNums((t-1)*m+1:t*m,:)*Ht(:,:,t)^(0.5);
    pseudorc(:,:,t) = r'*r;
end

%Truncate the data and the covariance
simulatedata=simulatedData(burnin+1:t,:);
pseudorc = pseudorc(:,:,burnin+1:t);
ht=Ht(:,:,burnin+1:t);

author: bashtage

repo: mfe-toolbox


text:
function ref = buildRefTraj(q,dq,ddq,ref)
% ref = buildRefTraj(q,dq,ddq,ref)
%
% This function computes the piecewise polynomial trajectories that make up
% the reference trajectory for the feedback linearization
%
% INPUTS:
%   q = [nConfig, nTime] = reference trajectory in configuration space
%   dq = [nConfig, nTime] = dq/dt = rates in time for ref traj
%   ddq = [nConfig, nTime] = ddq/ddt = accelerations for ref traj
%   ref = prototype for ref struct, with fields:
%       .wn = natural frequency of the controller
%       .xi = damping ratio of the controller
%       .c = [1, nConfig] = mapping from configuration to phase
%       .H = [nMeasure, nConfig] = mapping from configuration to measurement
%
% OUTPUTS:
%
%   ref = full ref struct, with added fields:
%       .pp = piecewise-polynomial trajectories for:
%           .h = reference measurement vector
%           .dh = dMeasurement/dPhase
%           .ddh = second derivative of measurements with respect to phase
%           .dhdt = dMeasurement/dTime
%
%

% Trajectories: measurement and phase vs time:
hRef = ref.H*q;   % Target measurement
dhRef = ref.H*dq;
ddhRef = ref.H*ddq;
pRef = ref.c*q; % Phase
dpRef = ref.c*dq;
ddpRef = ref.c*ddq;

% Compute derivatives wrt phase using chain rule:
dhRefdp = dhRef./dpRef;
ddhRefddp = (ddhRef - dhRefdp.*ddpRef)./(dpRef.^2);

% Represent trajectories as piecewise-polynomial
ref.pp.h = pchip(pRef,hRef);
ref.pp.dh = pchip(pRef,dhRefdp);   % Derivative wrt phase  (for ref traj)
ref.pp.dhdt = pchip(pRef,dhRef);   % Derivative wrt time   (for stabilization)
ref.pp.ddh = pchip(pRef,ddhRefddp);

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function CPD = set_clamped(CPD, bit)

CPD.clamped = bit;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function value = l4_to_i4 ( l )

%*****************************************************************************80
%
%% L4_TO_I4 converts an L4 to an I4.
%
%  Discussion:
%
%    0 is FALSE, and anything else if TRUE.
%
%    An I4 is an integer value.
%    An L4 is a logical value.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, logical L, a logical value.
%
%    Output, integer VALUE, the integer value of L.
%
  if ( l )
    value = 1;
  else
    value = 0;
  end

  return
end


text:
function xa = binom ( x, xx, npl, m, nt )

%*****************************************************************************80
%
%% BINOM: binomial expansion series for the (-1/M) power of a Chebyshev series.
%
%  Discussion:
%
%    This routine uses a certain number of terms of the binomial expansion
%    series to estimate the (-1/M) power of a given Chebyshev series.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Roger Broucke.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Roger Broucke,
%    Algorithm 446:
%    Ten Subroutines for the Manipulation of Chebyshev Series,
%    Communications of the ACM,
%    October 1973, Volume 16, Number 4, pages 254-256.
%
%  Parameters:
%
%    Input, real X(NPL), the given Chebyshev series.
%
%    Input, real XX(NPL), an initial estimate for
%    the Chebyshev series for the input function raised to the (-1/M) power.
%
%    Input, integer NPL, the number of terms in the
%    Chebyshev series.
%
%    Input, integer M, defines the exponent, (-1/M).
%    0 < M.
%
%    Input, integer NT, the number of terms of the binomial
%    series to be used.
%
%    Output, real XA(NPL), the estimated Chebyshev series
%    for the input function raised to the (-1/M) power.
%
  xa = zeros ( npl, 1 );

  dm = m;
  alfa = - 1.0 / dm;

  ww(1:npl) = x(1:npl);

  for k = 1 : m

    w2 = mltply ( ww, xx, npl );

    ww(1:npl) = w2(1:npl);

  end

  ww(1) = ww(1) - 2.0;

  xa(1) = 2.0;
  xa(2:npl) = 0.0;

  w3(1) = 2.0;
  w3(2:npl) = 0.0;

  for k = 2 : nt

    dkmm = k - 1;
    dkm2 = k - 2;
    coef = ( alfa - dkm2 ) / dkmm;

    w2 = mltply ( w3, ww, npl );

    w3(1:npl) = w2(1:npl) * coef;

    xa(1:npl) = xa(1:npl) + w3(1:npl);

  end

  w2 = mltply ( xa, xx, npl );

  xa(1:npl) = w2(1:npl);

  return
end


text:
function eeg_plot_metric

% eeg_plot_metric - create TeX plot label (Y is uV, X is msec)

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:52 $

% Licence:  GNU GPL, no express or implied warranties
% History:  07/2000, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Font.FontName   = 'Helvetica';
Font.FontUnits  = 'Pixels';
Font.FontSize   = 12;
Font.FontWeight = 'normal';
Font.FontAngle  = 'normal';

ylab = get(gca,'ylabel');
set(ylab,'string','\muV','Rotation',0,Font,'units','normalized',...
    'Position',[-0.10  0.95 0]);
xlab = get(gca,'xlabel');
set(xlab,'string','msec','Rotation',0,Font,'units','normalized',...
    'Position',[ 1.07 -0.01 0]);

return

author: PatternRecognition

repo: OpenBMI


text:
function phi = Segment2D_public( wb, wr, lambda, phi0, Nit )
% wb>0 edge detector
% wr region term
% lambda weight of region term over edge term
% phi0 initialization
% Nit number of iterations

verb = 1; fig_it = figure;
[Nx,Ny] = size(wb);
N = Nx*Ny;

% Parameters of augemented Lagrangian and minimizers
r1 = 1;
r2 = 1;
r3 = 0.1;
r4 = 0.1;
Njvarphi = 3; % increase it to get a better approximation to the minimizar with respect to varphi
Njphi = 2; % increase it to get better approximation to | grad phi | = 1
epsilon = 1; % controls smoothness of approximation to Heaviside and Dirac

% usufelu for FFT-based minimization
[Y,X] = meshgrid(0:Ny-1,0:Nx-1);
auxFFT = cos(2*pi/Nx*X)+cos(2*pi/Ny*Y)-2;

% initialize Lagrange multipliers and splitting variables
l1 = zeros(Nx,Ny);
l2x = zeros(Nx,Ny);
l2y = zeros(Nx,Ny);
l3 = zeros(Nx,Ny);
l4x = zeros(Nx,Ny);
l4y = zeros(Nx,Ny);
phi = phi0;
varphi = phi0;
u = H(phi0,epsilon);
dx = Fx(u);
dy = Fy(u);
px = Fx(phi);
py = Fy(phi);

wr = -lambda* wr;
% Start iterations: minimize w.r.t. each variable and update lagrange multipliers
for i=1:Nit

    % u^n+1 solved in Fourier domain
    b = r1*(H(varphi,epsilon)- l1/r1) - wr -r2* ( Bx(dx+l2x/r2)+By(dy+l2y/r2) );
    u = real(ifft2( fft2( b )./( -2*r2*auxFFT + r1 )));

    % d^n+1 solved by shrinkage
    Fxu = Fx(u);
    Fyu = Fy(u);
    [dx,dy] = shrink(Fxu - l2x/r2,Fyu - l2y/r2,wb/r2);

    % varphi^n+1 approx initialization and Newton method
    v = u + l1/r1;
    z = phi- l3/r3;
    E_0 = 0.5*r3*z.^2 + 0.5*r1*( 0.5 - v ).^2;
    E_z = 0.5*r1*( H(z,epsilon) - v ).^2;
        % initialize with solution to ideal Heaviside
    cond_m2 = ( z>=0 ).*( v>=0.5 ) + ( z<0 ).*( v<0.5 );
    varphi = z.*cond_m2 + (1-cond_m2).*( E_z < E_0 ).*z;
        % newton method to find zeros of derivative
    for j=1:Njvarphi
         df = r3*( varphi - z) + r1*( H(varphi,epsilon) - v ).*D(varphi,epsilon);
         hf = r3 + r1*D(varphi,epsilon).*D(varphi,epsilon) + r1*(H(varphi,epsilon) - v ).*dD(varphi,epsilon);
         varphi = varphi - df./( hf + (abs(hf)<1e-3) );
    end

    % lambda^n+1 update lagrangians
    l1 = l1 + r1* (u-H(varphi,epsilon));
    l2x = l2x + r2* (dx-Fxu);
    l2y = l2y + r2* (dy-Fyu);
    l3 = l3 + r3* (varphi- phi);

    %%%%%%%%%%%%%%%%%%%%%%%%%
    for jphi=1:Njphi

        % p^n+1 solved projecting in sphere
        ex = Fx(phi)-l4x/r4;
        ey = Fy(phi)-l4y/r4;
        norm_e = sqrt( ex.^2 + ey.^2 );
        px = ex./(norm_e + (norm_e==0));
        py = ey./(norm_e + (norm_e==0));

        % phi^n+1 solved in Fourier domain
        b = r3* (varphi+l3/r3) -r4.* ( Bx(px+l4x/r4)+By(py+l4y/r4) );
        phi = real(ifft2( fft2( b )./( -2*r4.*auxFFT + r3 )));

        % lambda^n+1 update lagrangians
        l4x = l4x + r4.*(px-Fx(phi));
        l4y = l4y + r4.*(py-Fy(phi));
    end

    %plotting
    if ((rem(i,5)==0)&&(verb==1))
        pause(0.01);
        figure(fig_it)
        clf; imagesc(phi); colormap(gray); colorbar;
        title([num2str(i), ' iterations']); hold on; contour(phi,[0 0],'m');
        for ct=-2:2; contour(phi,[ct ct],'k'); end; contour(phi,[0 0],'m');
    end
end

return;

% shrinkage function
function [xs,ys] = shrink(x,y,lambda)

s = sqrt(x.*conj(x)+y.*conj(y));
ss = s-lambda;
ss = ss.*(ss>0);

s = s+(s<lambda);
ss = ss./s;

xs = ss.*x;
ys = ss.*y;

return;

% smooth approximation of Heaviside function
function v = H(phi,epsilon)
v = 0.5*( 1 + 2/pi*atan(phi/epsilon) );
return;

% smooth approximation of Dirac distribution
function v = D(phi,epsilon)
v = 1/pi* epsilon./ (epsilon^2 + phi.^2);
return;

% smooth approximation of derivative of the Dirac distribution
function v = dD(phi,epsilon)
v = -2/pi* epsilon*phi./ (epsilon^2 + phi.^2).^2;

% Forward derivative operator on x with boundary condition u(:,1,:)=u(:,Nx,:)
function Fxu = Fx(u)
[Ny,Nx] = size(u);
Fxu = circshift(u,[0 -1])-u;
Fxu(:,Nx) = zeros(Ny,1);

% Forward derivative operator on y with boundary condition u(1,:,:)=u(Ny,:,:)
function Fyu = Fy(u)
[Ny,Nx] = size(u);
Fyu = circshift(u,[-1 0])-u;
Fyu(Ny,:) = zeros(1,Nx);

% Backward derivative operator on x with boundary condition Bxu(:,1)=u(:,1)
function Bxu = Bx(u)
[Ny,Nx] = size(u);
Bxu = u - circshift(u,[0 1]);
Bxu(:,1) = u(:,1);
Bxu(:,Nx) = -u(:,Nx-1);

% Backward derivative operator on x with boundary condition Bxu(1,:)=u(1,:)
function Byu = By(u)
[Ny,Nx] = size(u);
Byu = u - circshift(u,[1 0]);
Byu(1,:) = u(1,:);
Byu(Ny,:) = -u(Ny-1,:);


text:
function [mat_data, dims, folder_analysis] = distribute_data(nam, patch_dims, ...
    w_overlap, memory_size_per_patch, memory_size_to_use, dir_output, filter_kernel)
% divide the field of view into multiple small patches and save the data
% into these small patches individually. There are overlaps between these
% patches.

%% default values
if ~exist(nam, 'file')
    error('File does not exist!');
end
file_name = get_fullname(nam);

if ~exist('patch_dims', 'var')
    patch_dims = [];
end
if ~exist('w_overlap', 'var') || isempty(w_overlap)
    w_overlap = 10;
end
if ~exist('memory_size_per_patch', 'var') || isempty(memory_size_per_patch)
    memory_size_per_patch = 10;
end
if ~exist('memory_size_to_use', 'var') || isempty(memory_size_to_use)
    memory_size_to_use = 30;
end
if ~exist('filter_kernel', 'var') || isempty(filter_kernel)
    filter_kernel = [];
    filter_data = false;
else
    filter_data = true;
end

%% get the data dimension and determine the best way of distributing data
dims = get_data_dimension(file_name);
d1 = dims(1); d2 = dims(2); T = dims(3);

fprintf('\nThe data has %d X %d pixels X %d frames. \nLoading all data (double precision) requires %.3f GB RAM\n\n', d1, d2, T, prod(dims)/(2^27));

max_elements = memory_size_per_patch* (500^3); % x GB data can save x*1000^3/8 dobule numbers.
min_patch_width = 2*w_overlap+3;
max_patch_width = max(round(sqrt(double(round(max_elements/T))))-w_overlap*2, min_patch_width);

if isempty(patch_dims)
    % find the optimial batch size
    patch_dims = [1, 1] * patch_width;
else
    patch_dims(patch_dims< min_patch_width) = min_patch_width;
%     patch_dims(patch_dims>max_patch_width) = max_patch_width;
    if length(patch_dims)==1
        patch_dims = patch_dims * [1,1];
    else
        patch_dims = patch_dims(1:2);
    end
end

patch_dims = reshape(patch_dims, 1, []);
nr_patch = round(d1/patch_dims(1));
nc_patch = round(d2/patch_dims(2));
if nr_patch <= 1
    patch_idx_r = [1, d1];
else
    patch_idx_r = ceil(linspace(1, d1, nr_patch+1));
    patch_idx_r(end) = d1;
    if diff(patch_idx_r(1:2))<min_patch_width
        patch_idx_r = 1:min_patch_width:d1;
        patch_idx_r(end) = d1;
    end
end
nr_patch = length(patch_idx_r) - 1;

if nc_patch <= 1
    patch_idx_c = [1, d2];
else
    patch_idx_c = ceil(linspace(1, d2, nc_patch+1));
    if diff(patch_idx_c(1:2))<min_patch_width
        patch_idx_c = 1:min_patch_width:d2;
        patch_idx_c(end) = d2;
    end
end
nc_patch = length(patch_idx_c) - 1;

fprintf('The FOV is divided into %d X %d patches. \nEach patch has %d X %d pixels. \n',...
    nr_patch, nc_patch, diff(patch_idx_r(1:2)), diff(patch_idx_c(1:2)));
temp = prod(patch_dims+2*w_overlap)*T/(2^27);
fprintf('It requires %.3f GB RAM for loading data related to each patch. \n\n', temp);
if temp > memory_size_per_patch
    fprintf('You assigned a smaller memory for each patch.\n');
    fprintf('We suggest you process data in batch mode.\nEach batch has frames fewer than %d. \n\n', round(T*memory_size_per_patch/temp));
end

%% indices for dividing the FOV into multiple blocks.
block_idx_r = bsxfun(@plus, patch_idx_r, [-1-w_overlap; w_overlap]);
block_idx_c = bsxfun(@plus, patch_idx_c, [-1-w_overlap; w_overlap]);
block_idx_r = block_idx_r(:);
block_idx_c = block_idx_c(:);
block_idx_r(block_idx_r<1) = 1; block_idx_r(block_idx_r>d1) = d1;
block_idx_c(block_idx_c<1) = 1; block_idx_c(block_idx_c>d2) = d2;
block_idx_r = sort(unique(block_idx_r));
block_idx_c = sort(unique(block_idx_c));
nr_block = length(block_idx_r) - 1;
nc_block = length(block_idx_c) - 1;
% num_blocks = nr_block * nc_block;

%% prepare for distributing the data into multiple patches.
[dir_nm, file_nm, ~] = fileparts(file_name);
if ~exist('dir_output', 'var') || ~exist(dir_output, 'dir')
    dir_output = [dir_nm, filesep];
elseif dir_output(end)~=filesep
    dir_output = dir_output(1:(end-1));
end
folder_analysis = [dir_output, file_nm, '_source_extraction'];

mat_file = [folder_analysis, filesep,...
    sprintf('data_%d_%d_%d.mat', patch_dims(1), patch_dims(2), w_overlap)];

if ~exist(folder_analysis, 'dir')
    mkdir(folder_analysis);
end

% check whether the mat data has been created
if exist(mat_file, 'file')
    mat_data = matfile(mat_file);
    if all(patch_dims==mat_data.patch_dims) && (w_overlap==mat_data.w_overlap)
        fprintf('The data has been divided into multiple patches and saved into \n%s\n\n', mat_file);
        return;
    end
end

% create a mat file and save data info into it.
mat_data = matfile(mat_file, 'Writable', true);
save(mat_file, 'file_name', '-v7.3');
mat_data.patch_idx_r = patch_idx_r;
mat_data.patch_idx_c = patch_idx_c;
mat_data.nr_patch = nr_patch;
mat_data.nc_patch = nc_patch;
mat_data.block_idx_r = block_idx_r;
mat_data.block_idx_c = block_idx_c;
mat_data.nr_block = nr_block;
mat_data.nc_block = nc_block;
mat_data.w_overlap = w_overlap;
mat_data.patch_dims = patch_dims;
mat_data.dims = dims;

% pre-allocate matrices for saving the video data
img = smod_bigread2(file_name,1,1);
default_value = 0*img(1);
temp = whos('img');
mat_data.dtype = temp.class;
for m=1:nr_block
    r0 = block_idx_r(m);
    r1 = block_idx_r(m+1);
    nr = r1-r0+1;
    for n=1:nc_block
        c0 = block_idx_c(n);
        c1 = block_idx_c(n+1);
        nc = c1-c0+1;
        eval(sprintf('mat_data.Y_%d_%d_%d_%d(%d, %d, %d)=default_value;', r0, r1, c0, c1, nr, nc, T));
    end
end

% positions for determing a batch and the minimum block that includes the
% patch
patch_pos = cell(nr_patch, nc_patch);
block_pos = cell(nr_patch, nc_patch);
for m=1:nr_patch
    for n=1:nc_patch
        patch_pos{m, n} = [patch_idx_r(m:(m+1))+[0, -1*(m~=nr_patch)], patch_idx_c(n:(n+1))+[0, -1*(n~=nc_patch)]];
        block_pos{m, n} = [max(1, patch_idx_r(m)-w_overlap-1), min(d1, patch_idx_r(m+1)+w_overlap), ...
            max(1, patch_idx_c(n)-w_overlap-1), min(d2, patch_idx_c(n+1)+w_overlap)];
    end
end
mat_data.patch_pos = patch_pos;
mat_data.block_pos = block_pos;

%% load and distribute data
Tchunk = floor(memory_size_to_use * (500^3) / (d1*d2));
t_start= 0;
fprintf('\n-------- Loading --------\n');
fprintf('Data is being loaded and distributed into multiple small blocks for easy access.\n\n');

while t_start<T
    num2read = min(Tchunk, T-t_start);
    Y = smod_bigread2(file_name, t_start+1, num2read);
    Y(isnan(Y)) = 0;
    temp = Y(1);
    if filter_data
        Y = cast(imfilter(double(Y), filter_kernel, 'replicate'), 'like', temp);
    end
    for m=1:nr_block
        r0 = block_idx_r(m);
        r1 = block_idx_r(m+1);
        nr = r1-r0+1; %#ok<*NASGU>
        for n=1:nc_block
            c0 = block_idx_c(n);
            c1 = block_idx_c(n+1);
            nc = c1-c0+1;
            eval(sprintf('mat_data.Y_%d_%d_%d_%d(1:nr, 1:nc, %d:%d)=Y(%d:%d, %d:%d, :); ', ...
                r0, r1, c0, c1, t_start+1, t_start+num2read, r0, r1, c0, c1));
            fprintf('block(%2d,%2d)/(%d, %d) done\n', m, n, nr_block, nc_block);
        end
    end
    t_start = t_start + num2read;
    fprintf('loaded %d out of %d frames\n', t_start, T);
end
fprintf('\nThe data has been saved into \n%s\n', mat_file);
fprintf('\n-------- Done --------\n');

author: zhoupc


text:
% Author: Ricardo Baptista and Matthias Poloczek
% Date:   June 2018
%
% See LICENSE.md for copyright information
%

function [binary_sample, p] = logistic_sample(n_vars, logistic_A)
% LOGISTIC_SAMPLE: Function computes binary sample from the implemented
% logistic regression method, as well as the probability of the sample
% (likelihood under the binary logistic regression proposal).

% Declare new binary sample
binary_sample = zeros(1,n_vars);

% Declare new constant p
p = 1;

% Generate sample
for j=1:n_vars

    % Evaluate logistic regression function
    if j==1
        logistic_regressor = logistic_A(j,j);
    else
        logistic_regressor = logistic_A(j,j) + ...
                             logistic_A(j,1:j-1)*binary_sample(1:j-1)';
    end
    logistic_function = 1/(1 + exp(-logistic_regressor));

    % Generate sample from a bernoulli distribution
    binary_sample(j) = binornd(1,logistic_function);

    % Compute next success probability
    if (binary_sample(j) == 1)
        p = p*logistic_function;
    elseif (binary_sample(j) == 0)
        p = p*(1-logistic_function);
    end

end

end

author: baptistar

repo: BOCS

save_path: github-repos/MATLAB/baptistar-BOCS


text:
#!/usr/bin/env octave
%% Machine Learning Online Class
%  Exercise 7 | Principle Component Analysis and K-Means Clustering
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  exercise. You will need to complete the following functions:
%
%     pca.m
%     projectData.m
%     recoverData.m
%     computeCentroids.m
%     findClosestCentroids.m
%     kMeansInitCentroids.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% ================== Part 1: Load Example Dataset  ===================
%  We start this exercise by using a small dataset that is easily to
%  visualize
%
fprintf('Visualizing example dataset for PCA.\n\n');

%  The following command loads the dataset. You should now have the
%  variable X in your environment
load ('ex7data1.mat');

%  Visualize the example dataset
plot(X(:, 1), X(:, 2), 'bo');
axis([0.5 6.5 2 8]); axis square;

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =============== Part 2: Principal Component Analysis ===============
%  You should now implement PCA, a dimension reduction technique. You
%  should complete the code in pca.m
%
fprintf('\nRunning PCA on example dataset.\n\n');

%  Before running PCA, it is important to first normalize X
[X_norm, mu, sigma] = featureNormalize(X);

%  Run PCA
[U, S] = pca(X_norm);

%  Compute mu, the mean of the each feature

%  Draw the eigenvectors centered at mean of data. These lines show the
%  directions of maximum variations in the dataset.
hold on;
drawLine(mu, mu + 1.5 * S(1,1) * U(:,1)', '-k', 'LineWidth', 2);
drawLine(mu, mu + 1.5 * S(2,2) * U(:,2)', '-k', 'LineWidth', 2);
hold off;

fprintf('Top eigenvector: \n');
fprintf(' U(:,1) = %f %f \n', U(1,1), U(2,1));
fprintf('\n(you should expect to see -0.707107 -0.707107)\n');

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =================== Part 3: Dimension Reduction ===================
%  You should now implement the projection step to map the data onto the
%  first k eigenvectors. The code will then plot the data in this reduced
%  dimensional space.  This will show you what the data looks like when
%  using only the corresponding eigenvectors to reconstruct it.
%
%  You should complete the code in projectData.m
%
fprintf('\nDimension reduction on example dataset.\n\n');

%  Plot the normalized dataset (returned from pca)
plot(X_norm(:, 1), X_norm(:, 2), 'bo');
axis([-4 3 -4 3]); axis square

%  Project the data onto K = 1 dimension
K = 1;
Z = projectData(X_norm, U, K);
fprintf('Projection of the first example: %f\n', Z(1));
fprintf('\n(this value should be about 1.481274)\n\n');

X_rec  = recoverData(Z, U, K);
fprintf('Approximation of the first example: %f %f\n', X_rec(1, 1), X_rec(1, 2));
fprintf('\n(this value should be about  -1.047419 -1.047419)\n\n');

%  Draw lines connecting the projected points to the original points
hold on;
plot(X_rec(:, 1), X_rec(:, 2), 'ro');
for i = 1:size(X_norm, 1)
    drawLine(X_norm(i,:), X_rec(i,:), '--k', 'LineWidth', 1);
end
hold off

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =============== Part 4: Loading and Visualizing Face Data =============
%  We start the exercise by first loading and visualizing the dataset.
%  The following code will load the dataset into your environment
%
fprintf('\nLoading face dataset.\n\n');

%  Load Face dataset
load ('ex7faces.mat')

%  Display the first 100 faces in the dataset
displayData(X(1:100, :));

fprintf('Program paused. Press enter to continue.\n');
pause;

%% =========== Part 5: PCA on Face Data: Eigenfaces  ===================
%  Run PCA and visualize the eigenvectors which are in this case eigenfaces
%  We display the first 36 eigenfaces.
%
fprintf(['\nRunning PCA on face dataset.\n' ...
         '(this mght take a minute or two ...)\n\n']);

%  Before running PCA, it is important to first normalize X by subtracting
%  the mean value from each feature
[X_norm, mu, sigma] = featureNormalize(X);

%  Run PCA
[U, S] = pca(X_norm);

%  Visualize the top 36 eigenvectors found
displayData(U(:, 1:36)');

fprintf('Program paused. Press enter to continue.\n');
pause;

%% ============= Part 6: Dimension Reduction for Faces =================
%  Project images to the eigen space using the top k eigenvectors
%  If you are applying a machine learning algorithm
fprintf('\nDimension reduction for face dataset.\n\n');

K = 100;
Z = projectData(X_norm, U, K);

fprintf('The projected data Z has a size of: ')
fprintf('%d ', size(Z));

fprintf('\n\nProgram paused. Press enter to continue.\n');
pause;

%% ==== Part 7: Visualization of Faces after PCA Dimension Reduction ====
%  Project images to the eigen space using the top K eigen vectors and
%  visualize only using those K dimensions
%  Compare to the original input, which is also displayed

fprintf('\nVisualizing the projected (reduced dimension) faces.\n\n');

K = 100;
X_rec  = recoverData(Z, U, K);

% Display normalized data
subplot(1, 2, 1);
displayData(X_norm(1:100,:));
title('Original faces');
axis square;

% Display reconstructed data from only k eigenfaces
subplot(1, 2, 2);
displayData(X_rec(1:100,:));
title('Recovered faces');
axis square;

fprintf('Program paused. Press enter to continue.\n');
pause;

%% === Part 8(a): Optional (ungraded) Exercise: PCA for Visualization ===
%  One useful application of PCA is to use it to visualize high-dimensional
%  data. In the last K-Means exercise you ran K-Means on 3-dimensional
%  pixel colors of an image. We first visualize this output in 3D, and then
%  apply PCA to obtain a visualization in 2D.

close all; close all; clc

% Re-load the image from the previous exercise and run K-Means on it
% For this to work, you need to complete the K-Means assignment first
A = double(imread('bird_small.png'));

% If imread does not work for you, you can try instead
%   load ('bird_small.mat');

A = A / 255;
img_size = size(A);
X = reshape(A, img_size(1) * img_size(2), 3);
K = 16;
max_iters = 10;
initial_centroids = kMeansInitCentroids(X, K);
[centroids, idx] = runkMeans(X, initial_centroids, max_iters);

%  Sample 1000 random indexes (since working with all the data is
%  too expensive. If you have a fast computer, you may increase this.
sel = floor(rand(1000, 1) * size(X, 1)) + 1;

%  Setup Color Palette
palette = hsv(K);
colors = palette(idx(sel), :);

%  Visualize the data and centroid memberships in 3D
figure;
scatter3(X(sel, 1), X(sel, 2), X(sel, 3), 10, colors);
title('Pixel dataset plotted in 3D. Color shows centroid memberships');
fprintf('Program paused. Press enter to continue.\n');
pause;

%% === Part 8(b): Optional (ungraded) Exercise: PCA for Visualization ===
% Use PCA to project this cloud to 2D for visualization

% Subtract the mean to use PCA
[X_norm, mu, sigma] = featureNormalize(X);

% PCA and project the data to 2D
[U, S] = pca(X_norm);
Z = projectData(X_norm, U, 2);

% Plot in 2D
figure;
plotDataPoints(Z(sel, :), idx(sel), K);
title('Pixel dataset plotted in 2D, using PCA for dimensionality reduction');
fprintf('Program paused. Press enter to continue.\n');
pause;

author: SaveTheRbtz

repo: ml-class


text:
function [pop] = initializega(num, bounds, evalFN,evalOps,options)
% function [pop]=initializega(populationSize, variableBounds,evalFN,
%                           evalOps,options)
%    initializega creates a matrix of random numbers with
%    a number of rows equal to the populationSize and a number
%    columns equal to the number of rows in bounds plus 1 for
%    the f(x) value which is found by applying the evalFN.
%    This is used by the ga to create the population if it
%    is not supplied.
%
% pop            - the initial, evaluated, random population
% populatoinSize - the size of the population, i.e. the number to create
% variableBounds - a matrix which contains the bounds of each variable, i.e.
%                  [var1_high var1_low; var2_high var2_low; ....]
% evalFN         - the evaluation fn, usually the name of the .m file for
%                  evaluation
% evalOps        - any options to be passed to the eval function defaults []
% options        - options to the initialize function, ie.
%                  [type prec] where eps is the epsilon value
%                  and the second option is 1 for float and 0 for binary,
%                  prec is the precision of the variables defaults [1e-6 1]

% Binary and Real-Valued Simulation Evolution for Matlab GAOT V2
% Copyright (C) 1998 C.R. Houck, J.A. Joines, M.G. Kay
%
% C.R. Houck, J.Joines, and M.Kay. A genetic algorithm for function
% optimization: A Matlab implementation. ACM Transactions on Mathmatical
% Software, Submitted 1996.
%
% 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 1, 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. A copy of the GNU
% General Public License can be obtained from the
% Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

if nargin<5
  options=[1e-6 1];
end
if nargin<4
  evalOps=[];
end

if any(evalFN<48) %Not a .m file
  if options(2)==1 %Float GA
    estr=['x=pop(i,1); pop(i,xZomeLength)=', evalFN ';'];
  else %Binary GA
    estr=['x=b2f(pop(i,:),bounds,bits); pop(i,xZomeLength)=', evalFN ';'];
  end
else %A .m file
  if options(2)==1 %Float GA
    estr=['[ pop(i,:) pop(i,xZomeLength)]=' evalFN '(pop(i,:),[0 evalOps]);'];
  else %Binary GA
    estr=['x=b2f(pop(i,:),bounds,bits);[x v]=' evalFN ...
	'(x,[0 evalOps]); pop(i,:)=[f2b(x,bounds,bits) v];'];
    end
end

numVars     = size(bounds,1); 		%Number of variables
rng         = (bounds(:,2)-bounds(:,1))'; %The variable ranges'

if options(2)==1 %Float GA
  xZomeLength = numVars+1; 		%Length of string is numVar + fit
  pop         = zeros(num,xZomeLength); 	%Allocate the new population
  pop(:,1:numVars)=(ones(num,1)*rng).*(rand(num,numVars))+...
    (ones(num,1)*bounds(:,1)');
else %Binary GA
  bits=calcbits(bounds,options(1));
  xZomeLength = sum(bits)+1; 		%Length of string is numVar + fit
  pop = round(rand(num,sum(bits)+1));
end

for i=1:num
  eval(estr);
end

author: HuangCongQing


text:
% Calculates DNA-bound protein displacements
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 4/6/2012
function runDisplacementCalculation(simBatch, simIdx, outDir)
setWarnings;
setPath;

%import classes
import edu.stanford.covert.cell.sim.analysis.Figures34;
import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;

%parse inputs
if ischar(simIdx)
	simIdx = str2num(simIdx);
end

%run calculation
simBatchDir = [SimulationDiskUtil.getBaseDir() filesep simBatch filesep];
Figures34.calculateDNABoundProteinDisplacements(outDir, simBatchDir, simIdx, 1, 1);

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [dynpcm2] = atm2dynpcm2(atm)
% Convert pressure from atmospheres to dynes per square centimeter.
%
% So I watched that movie "Winnebago Man" last night.  It was pretty good.
% I didn't expect much from a low-budget documentary about the guy from
% some novelty viral internet videos to be very good or entertaining.  And
% I expected the first half the be good, then the second half to drag, as
% happens to so many indie docs.  But in fact it was funny at the
% beginning, then became MORE engaging in the second half.  That was a nice
% surprise.
% Chad A. Greene
dynpcm2 = atm*1.01325e+6;

author: Sable


text:
function output = F_inner_product(input_layers)

input1 = input_layers{1}.a;
input2 = input_layers{2}.a;

output = sum(input1 .* input2);

end

author: singaxiong

repo: SignalGraph


text:
function element_size = p04_element_size ( )

%*****************************************************************************80
%
%% P04_ELEMENT_SIZE returns a typical element size for problem 04.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, real ELEMENT_SIZE, a typical element size.
%
  element_size = 0.1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function ival = c8mat_is_unitary ( m, n, a )

%*****************************************************************************80
%
%% C8MAT_IS_UNITARY checks whether a complex matrix is unitary.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the row and column dimensions
%    of the matrix.  M and N must be positive.
%
%    Input, complex A(M,N), the matrix.
%
%    Output, integer IVAL:
%    -1, the matrix is not unitary; M /= N.
%    -2, the matrix is not unitary; row I * column I /= 1;
%    -3, the matrix is not unitary; row I * column J /= 0;
%    1, the matrix is unitary.
%
  tol = 0.0001;

  if ( m ~= n )
    ival = -1;
    return
  end

  for i = 1 : n

    test = a(i,1:n) * conj ( a(1:n,i) ) - 1.0;

    if ( tol < abs ( test ) )
      ival = -2;
      return
    end

    for j = i + 1 : n

      test = a(i,1:n ) * conj ( a(1:n,j) );

      if ( tol < abs ( test ) )
        ival = -3;
        return
      end

    end

  end

  ival = 1;

  return
end


text:
classdef RunningVademecumInParalel < handle

    properties (Access = private)
        samplePoints
        cellVar
        fileName
    end

    properties (Access = private)
        nCase
        nMx
        nMy
        nPhi
        phiMin
        phiMax
    end

    methods (Access = public)

        function obj = RunningVademecumInParalel(cParams)
            obj.init(cParams);
            obj.compute();
        end

    end

     methods (Access = private)

        function init(obj,nC)
            obj.nCase  = nC;
            obj.nMx    = 20;
            obj.nMy    = 20;
            obj.nPhi   = 10;
        end

        function compute(obj)
            obj.computeSameStressSign();
            obj.computeDifferentStressSign();
        end

        function computeSameStressSign(obj)
            obj.phiMin = 0;
            obj.phiMax = pi/4;
            obj.fileName = 'OptimalSuperEllipseSameStressSign';
            obj.computeVademecums();
        end

        function computeDifferentStressSign(obj)
            obj.phiMin = pi/2;
            obj.phiMax = 3*pi/4;
            obj.fileName = 'OptimalSuperEllipseDifferentStressSign';
            obj.computeVademecums();
        end

        function computeVademecums(obj)
            obj.createSamplePoints();
            obj.cellVar = cell(obj.nMx,obj.nMy,obj.nPhi);
            imx = obj.nCase;
            for imy = 1:obj.nMy
                obj.displayPercentatge(imx,imy);
                for iphi = 1:obj.nPhi
                    c = obj.computeOneVademecum(imx,imy,iphi);
                    obj.cellVar{imx,imy,iphi} = c;
                end
            end
            c = obj.cellVar;
            save([obj.fileName,num2str(imx)],'c');
        end

        function cellV = computeOneVademecum(obj,imx,imy,iphi)
            iter = obj.computeGlobalIteration(imx,imy);
            s.rho = obj.samplePoints.rhoV(iter);
            s.xi = obj.samplePoints.txiV(iter);
            s.phi = obj.samplePoints.phiV(iphi);
            s.nCase = [num2str(imx),num2str(imy),num2str(iphi)];
            v = VademecumCalculator(s);
            v.compute();
            %v.cellVariables = s.nCase;
            cellV = v.cellVariables;
        end

       function createSamplePoints(obj)
            s.type = 'FromMxMy';
            s.nMx  = obj.nMx;
            s.nMy  = obj.nMy;
            s.nPhi = obj.nPhi;
            s.phiMin = obj.phiMin;
            s.phiMax = obj.phiMax;
            sample = SamplePointsCreatorForOptimalExponentComputer.create(s);
            sample.compute();
            obj.samplePoints = sample;
        end

        function iter = computeGlobalIteration(obj,imx,imy)
            iter = (imy + obj.nMx*(imx -1));
        end

        function displayPercentatge(obj,imx,imy)
            perc = imy/(obj.nMy)*100;
            %disp([num2str(perc),'% done']);
            fid = fopen('Output.txt','a+');
            fprintf(fid,'%6.2f %%\n',perc);
            fclose(fid);
        end

     end

end


text:
function variance = uniform_variance ( a, b )

%*****************************************************************************80
%
%% UNIFORM_VARIANCE returns the variance of the Uniform PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    A < B.
%
%    Output, real VARIANCE, the variance of the PDF.
%
  variance = ( b - a )^2 / 12.0;

  return
end


text:
n = 10;
N = 2^n;
K = 2^(n - 1);
R = K/N;
design_snr = 2.5;%dB
%You;d better to learn more about the relationship between RM code and
%polar code. Then you can construct good polar code.
sigma = 1/sqrt(2 * R) * 10^(-design_snr/20);
max_runs = 1e5;
llr_layer_vec = get_llr_layer(N);
bit_layer_vec = get_bit_layer(N);
lambda_offset = 2.^(0 : n);
ber = zeros(N, 1);
for i_run = 1 : max_runs
    if mod(i_run, max_runs/100) == 1
        disp(['Type I monte carlo code construction running = ' num2str(i_run/max_runs*100) '%']);
    end
    dummy_info = rand(N, 1) > 0.5;
    noise = randn(N, 1);
    x = my_polar_encode(dummy_info, lambda_offset, llr_layer_vec);
    bpsk = 1 - 2 * x;
    y = bpsk + sigma * noise;
    llr = 2 * y / sigma^2;
    ber_tmp = mc_typeI_SC_decoder(llr, lambda_offset, llr_layer_vec, bit_layer_vec, dummy_info);
    ber = ber + ber_tmp;
end
[~, channel_ordered] = sort(ber);
info_bits = sort(channel_ordered(1 : K));
disp('Type I MC CC done.')
disp('Variable "info_bits" is what you want')

author: YuYongRun

repo: PolarCodeDecodersInMatlab


text:
function fulpts = fancySelect(curImage,curSize,curMin, curMax, X, Y, zoomSize, selOff)
%
% pts = selectPoint
%    Select a set of points from the current image by clicking.
%    left mouse button -- selects a point
%    middle mouse button -- deletes a point
%    right mouse button -- exits
%
%    Started by B. Wandell, 06.28.93
%
% Variable Declarations
num = 0;		% The number of points selected
ord = zeros(curSize(1),curSize(2));
fulpts = zeros(curSize(1),curSize(2));

oldImage = curImage;		% For erasing
num = 1;
while 1 == 1
 [x y but] =  mrGinput(1,'cross');
  x = floor(x); y = floor(y);
  fulx = x + selOff(1) - 1;	% "-1" offset is required becuase the upper
  fuly = y + selOff(2) - 1;	% left hand corner of the screen is (1,1) not (0,0)
  if x <= zoomSize(2) & y <= zoomSize(1) & x > 0 & y > 0
     thept = fuly+(fulx-1)*curSize(1);
     if but == 1
       fulpts(fuly,fulx) = 1;
       ord(fuly,fulx) = num;
       curImage(thept) = -2;		% blue
       num = num + 1;
     elseif but == 2
       fulpts(fuly,fulx) = 0;
       curImage(thept) = oldImage(thept);
     elseif but == 3
       break
     end
  myShowImageVol(curImage, curSize, curMin, curMax, X, Y, zoomSize, selOff);
  else
   s = sprintf('Out of Range: %d %d',x,y)
  end
end

fulpts = reshape(fulpts,1,prod(size(fulpts)));
ord = reshape(ord,1,prod(size(fulpts))) .* fulpts;
dum = (ord == 0);
ord(dum) = 99999.*ones(sum(dum),1);
[srt, arang] = sort(ord);
fulpts = arang(1:sum(fulpts));

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ctheta_f, theta_pre_f, L] = train_filter(ctheta_f, xlf, yf, ...
    theta_pre_f, params, output_sz, seq)

    for k = 1: numel(xlf)
        x_f = xlf{k};

        if (seq.frame == 1)
            theta_pre_f{k} = zeros(size(x_f));
            lambda_2 = 0;
        else
            lambda_2 = params.lambda2(k);
        end

        % intialisation
        theta_f = single(zeros(size(x_f)));
        theta_prime_f = theta_f;
        eta_f = theta_f;
        mu  = params.init_penalty_factor(k);
        mu_scale_step = params.penalty_scale_step(k);

        % pre-compute the variables
        T = prod(output_sz);
        S_xx = sum(conj(x_f) .* x_f, 3);
        Stheta_pre_f = sum(conj(x_f) .* theta_pre_f{k}, 3);
        Sfx_pre_f = bsxfun(@times, x_f, Stheta_pre_f);

        % solve via ADMM algorithm
        iter = 1;
        while (iter <= params.max_iterations)

            % solving theta
            B = S_xx + T * (mu + lambda_2);
            Sgx_f = sum(conj(x_f) .* theta_prime_f, 3);
            Shx_f = sum(conj(x_f) .* eta_f, 3);

            theta_f = ((1/(T*(mu + lambda_2)) * bsxfun(@times,  yf{k}, x_f)) ...
                - ((1/(mu + lambda_2)) * eta_f) +(mu/(mu + lambda_2)) ...
                * theta_prime_f) + (lambda_2/(mu + lambda_2)) * theta_pre_f{k} ...
                - bsxfun(@rdivide,(1/(T*(mu + lambda_2)) * bsxfun(@times, ...
                x_f, (S_xx .*  yf{k})) + (lambda_2/(mu + lambda_2)) * Sfx_pre_f - ...
                (1/(mu + lambda_2))* (bsxfun(@times, x_f, Shx_f)) +(mu/(mu ...
                + lambda_2))* (bsxfun(@times, x_f, Sgx_f))), B);

            % solving theta_prime
            X = real(ifft2(mu * theta_f+ eta_f));
            if (seq.frame == 1)
                X_temp = zeros(size(X));
                for i = 1:size(X,3)
                 X_temp(:,:,i) = X(:,:,i) ./  (params.reg_window{k} .^2 + mu);
                end
                L = 0;
            else
            X_temp=X;
            L{k} = max(0,1-1./(mu*numel(X)*sqrt(sum(X_temp.^2,3))));

            [~,b] = sort(L{k}(:),'descend');
            L{k}(b(ceil(params.fs_rate(k)*1/params.search_area_scale^2*numel(b)):end)) = 0;

            X_temp = repmat(L{k},1,1,size(X_temp,3)) .* X_temp;
            end

            theta_prime_f = fft2(X_temp);

            %   update eta
            eta_f = eta_f + (mu * (theta_f - theta_prime_f));

            %   update mu
            mu = min(mu_scale_step * mu, 0.1);

            iter = iter+1;
        end

        % save the trained filters
        theta_pre_f{k} = theta_f;

        if seq.frame == 1
            ctheta_f{k} = theta_f;
        else
            ctheta_f{k} = params.rl * theta_f + (1-params.rl) * ctheta_f{k};
        end
    end

end

author: XU-TIANYANG

repo: LADCF

save_path: github-repos/MATLAB/XU-TIANYANG-LADCF


text:
function mono_total_next_grevlex_test ( )

%*****************************************************************************80
%
%% MONO_TOTAL_NEXT_GREVLEX_TEST tests MONO_TOTAL_NEXT_GREVLEX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 December 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MONO_TOTAL_NEXT_GREVLEX_TEST\n' );
  fprintf ( 1, '  MONO_TOTAL_NEXT_GREVLEX can list the monomials\n' );
  fprintf ( 1, '  in M variables, of total degree N,\n' );
  fprintf ( 1, '  one at a time, in graded reverse lexicographic order.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We start the process with (0,0,...,0,N).\n' );
  fprintf ( 1, '  The process ends with (N,0,...,0,0)\n' );

  n = 3;
  m = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Let M = %d\n', m );
  fprintf ( 1, '      N = %d\n', n );
  fprintf ( 1, '\n' );

  x = [ 0, 0, n ];
  i = 1;

  while ( 1 )

    fprintf ( 1, '  %2d:', i );
    for j = 1 : m
      fprintf ( 1, '  %1d', x(j) );
    end
    fprintf ( 1, '\n' );

    if ( x(1) == n )
      break
    end

    x = mono_total_next_grevlex ( m, n, x );
    i = i + 1;

  end

  return
end


text:
% PURPOSE  : Assign events to bins
%
% FORMAT   :
%
% EEG  = pop_binlister( EEG , parameters);
%
% INPUTS   :
%
% EEG           - input dataset
%
% The available parameters are as follows:
%
%     'BDF'         - name of the text file containing your bin descriptions (formulas).
%     'ImportEL'	  - (optional) name of the text file, to import, that contain the event information to process,
%                     according to ERPLAB format (see tutorial).
%
%     'ExportEL' 	  - (optional) name of the text file, to export, that will contain the upgraded event information,
%                     according to ERPLAB format (see tutorial).

%     'Resetflag'   - set (all) flags to zero before starting binlister process. 'on'=reset;  'off':keep as it is.
%
%     'Forbidden'	  - array of event codes (numeric). If any of these codes is among a set of codes successfully captured by a bin
%                     this "capture" will be disable.
%     'Ignore'      - array of event codes (numeric) to be ignored. Binlister will be blind to them.
%
%     'UpdateEEG'   - after binlister process you can move the upgraded event information to EEG.event field. 'on'=update, 'off'=keep as it is.
%     'Warning'     - 'on'- warn if EVENTLIST will be overwritten. 'off' - do not warn if EVENTLIST will be overwritten.
%     'SendEL2'     - once binlister ends its work, you can send a copy of the resulting EVENTLIST structure to:
%                    'Text'           - send to text file
%                    'EEG'            - send to EEG structure
%                    'EEG&Text'       - send to EEG & text file
%                    'Workspace'      - send to Matlab workspace,
%                    'Workspace&Text' - send to Workspace and text file,
%                    'Workspace&EEG'  - send to workspace and EEG,
%                    'All'- send to all of them.
%     'Report'      - 'on'= create report about binlister performance, 'off'= do not create a report.
%     'Saveas'      - (optional) open GUI for saving dataset. 'on'/'off'
%
%
%
% OUTPUTS  :
%
% EEG            - updated output dataset
%
% EXAMPLE  :
%
% EEG  = pop_binlister( EEG , 'BDF', '/Users/etfoo/Documents/MATLAB/Test_Data/binlister_demo_2.txt', 'ExportEL', ...
%                      '/Users/etfoo/Documents/MATLAB/Test_output.txt', 'Forbidden',  12, 'Ignore',  10, 'ImportEL',...
%                      '/Users/etfoo/Documents/MATLAB/Test.txt', 'Saveas','on', 'SendEL2', 'All', 'UpdateEEG', 'on', 'Warning', 'on');
%
%
% See also menuBinListGUI.m binlister.m
%
% *** This function is part of ERPLAB Toolbox ***
% Author: Javier Lopez-Calderon and Steven Luck
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% 2009

%b8d3721ed219e65100184c6b95db209bb8d3721ed219e65100184c6b95db209b
%
% ERPLAB Toolbox
% Copyright ? 2007 The Regents of the University of California
% Created by Javier Lopez-Calderon and Steven Luck
% Center for Mind and Brain, University of California, Davis,
% javlopez@ucdavis.edu, sjluck@ucdavis.edu
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 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/>.

function [EEG, com] = pop_binlister_fast(EEG, varargin)
global ALLEEG
global CURRENTSET
com = '';
if nargin < 1
        help pop_binlister
        return
end
if isobject(EEG) % eegobj
        whenEEGisanObject % calls a script for showing an error window
        return
end
if nargin==1
        if iseegstruct(EEG) && length(EEG)>1
                msgboxText = 'Unfortunately, binlister does not work with multiple datasets';
                title      = 'ERPLAB: multiple inputs';
                errorfound(msgboxText, title);
                return
        end
        if isfield(EEG.event, 'type')
                if ~all(cellfun(@isnumeric, { EEG.event.type }))
                        msgboxText = ['Some or all of your events contain only a text-based event label, '...
                                'and not a numeric event code.  ERPLAB must have a numeric code '...
                                'for every event (a text label can also be present).  For labels such as "S14", '...
                                'you can automatically create an equivalent numeric event code (e.g., 14) by '...
                                'checking the box labeled "Create numeric equivalents of nonnumeric event codes '...
                                'when possible" when creating the EventList.\n\nFor labels that cannot be automatically '...
                                'converted (e.g., "RESP"), or for ambiguous cases (e.g., when you have both "S14" and '...
                                '"R14"), you can flexibly create numeric versions using the Advanced button at EVENTLIST GUI.'];
                        title      = 'ERPLAB: Warning';
                        errorfound(sprintf(msgboxText), title);
                        return
                end
        end
        def  = erpworkingmemory('pop_binlister');
        if isempty(def)
                def = {'' '' '' 0 [] [] 0 0 0 1 0};
        end
        try
                ERPX = evalin('base', 'ERP');
        catch
                ERPX = [];
        end

        %
        % Call a GUI
        %
        packarray = menuBinListGUI(EEG, ERPX, def);

        if isempty(packarray)
                disp('User selected Cancel')
                return
        end

        file1      = packarray{1};         % bin descriptor file
        file2      = packarray{2};         % external eventlist (read event list from)
        file3      = packarray{3};         % text file containing the updated EVENTLIST (Write resulting eventlist to)
        flagrst    = packarray{4};         % 1 means reset flags
        forbiddenCodeArray = packarray{5};
        ignoreCodeArray    = packarray{6};
        updevent   = packarray{7};
        option2do  = packarray{8};         % See  option2do below
        reportable = packarray{9};         % 1 means create a report about binlister work.
        iswarning  = packarray{10};        % 1 means create a report about binlister work.
        getfromerp = packarray{11};
        indexEL    = packarray{12};

        if isempty(file2) || strcmpi(file2,'no') || strcmpi(file2,'none')
                if getfromerp
                        if isempty(ERPX)
                                msgboxText =  'You do not have any loaded ERPset yet.';
                                title = 'ERPLAB: no data';
                                errorfound(msgboxText, title);
                                return
                        end
                        EEGaux = EEG;
                        EEG    = ERPX;
                end
                if isfield(EEG, 'EVENTLIST')
                        %if length(EEG.EVENTLIST)>1
                        %        prompt    = {'Enter  EVENLIST index:'};
                        %        dlg_title = 'Multiple EVENTLIST at EEG were detected';
                        %        num_lines = 1;
                        %        def       = {'1'};
                        %        answer    = inputvalue(prompt,dlg_title,num_lines,def);
                        %        if isempty(answer)
                        %                disp('User selected Cancel')
                        %                return
                        %        end
                        %        indexEL = str2num(answer{1});
                        %else
                        %        indexEL = 1;
                        %end

                        if isfield(EEG.EVENTLIST, 'eventinfo')
                                if isempty(EEG.EVENTLIST(indexEL).eventinfo)
                                        msgboxText = ['EVENTLIST.eventinfo structure is empty!\n'...
                                                'Use Create EVENTLIST before BINLISTER'];
                                        title = 'ERPLAB: Error';
                                        errorfound(sprintf(msgboxText), title);
                                        return
                                end
                        else
                                msgboxText =  ['EVENTLIST.eventinfo structure was not found!\n'...
                                        'Use Create EVENTLIST before BINLISTER'];
                                title = 'ERPLAB: Error';
                                errorfound(sprintf(msgboxText), title);
                                return
                        end
                else
                        msgboxText =  ['EVENTLIST structure was not found!\n'...
                                'Use Create EVENTLIST before BINLISTER'];
                        title = 'ERPLAB: Error';
                        errorfound(sprintf(msgboxText), title);
                        return
                end
                logfilename = 'no';
                logpathname = '';
                file2 = [logpathname logfilename];
                disp('For LOGFILE, user selected INTERNAL')
        end

        erpworkingmemory('pop_binlister', {file1, file2, file3, flagrst, forbiddenCodeArray, ignoreCodeArray,...
                updevent, option2do, reportable, iswarning, getfromerp, indexEL});

        if flagrst==1
                strflagrst = 'on';
        else
                strflagrst = 'off';
        end
        if updevent==1
                strupdevent = 'on';
        else
                strupdevent = 'off';
        end
        switch option2do
                case 0
                        msgboxText = 'Where should I send the update EVENTLIST???\n Pick an option.';
                        title = 'ERPLAB: Error';
                        errorfound(sprintf(msgboxText), title);
                        return
                case 1
                        stroption2do = 'Text';
                case 2
                        stroption2do = 'EEG';
                case 3
                        stroption2do = 'EEG&Text';
                case 4
                        stroption2do = 'Workspace';
                case 5
                        stroption2do = 'Workspace&Text';
                case 6
                        stroption2do = 'Workspace&EEG';
                case 7
                        stroption2do = 'All';
        end
        if reportable==1
                strreportable = 'on';
        else
                strreportable = 'off';
        end
        if iswarning==1
                striswarning = 'on';
        else
                striswarning = 'off';
        end
        if ~getfromerp
                EEG.setname = [EEG.setname '_bins']; %suggest a new name
        end

        %
        % Somersault
        %
        [EEG, com] = pop_binlister(EEG, 'BDF', file1, 'ImportEL', file2, 'ExportEL', file3, 'Resetflag', strflagrst, 'Forbidden', forbiddenCodeArray,...
                'Ignore', ignoreCodeArray, 'UpdateEEG', strupdevent, 'SendEL2', stroption2do,'Report', strreportable, 'Warning', striswarning,...
                'Saveas', 'on', 'IndexEL', indexEL, 'History', 'gui');
        if getfromerp
                EEG = EEGaux;
        end
        return
end

%
% Parsing inputs
%
p = inputParser;
p.FunctionName  = mfilename;
p.CaseSensitive = false;
p.addRequired('EEG');
% option(s)
p.addParamValue('BDF', '', @ischar);
p.addParamValue('ImportEL', 'none', @ischar);
p.addParamValue('ExportEL', 'none', @ischar);
p.addParamValue('Resetflag', 'off', @ischar);% 'on', 'off'
p.addParamValue('Forbidden', [], @isnumeric);
p.addParamValue('Ignore', [], @isnumeric);
p.addParamValue('UpdateEEG', 'off', @ischar);% 'on', 'off'
p.addParamValue('SendEL2', 'EEG', @ischar);  %
p.addParamValue('Report', 'off', @ischar);   % 'on', 'off'
p.addParamValue('Warning', 'off', @ischar);  % 'on', 'off'
p.addParamValue('Saveas', 'off', @ischar);   % 'on', 'off'
p.addParamValue('IndexEL', 1, @isnumeric);   % integer
p.addParamValue('Voutput', 'EEG', @ischar);  % 'EEG' (which contains EVENTLIST) or 'EVENTLIST' (EVENTLIST only)
p.addParamValue('History', 'script', @ischar); % history from scripting

p.parse(EEG, varargin{:});

if length(EEG)>1
        msgboxText =  'ERPLAB says: Unfortunately, this function does not work with multiple datasets';
        error(msgboxText);
end

file1      = p.Results.BDF;              % bin descriptor file
file2      = p.Results.ImportEL;         % external eventlist (read event list from)
file3      = p.Results.ExportEL;         % text file containing the updated EVENTLIST (Write resulting eventlist to)
indexEL    = p.Results.IndexEL;          % EVENTLIST's index (in case of multiple EVENTLISTs)
outputv    = p.Results.Voutput;          % 'EEG' (which contains EVENTLIST) or 'EVENTLIST' (EVENTLIST only)

if strcmpi(p.Results.Resetflag, 'on')
        flagrst = 1;
else
        flagrst = 0;
end
if strcmpi(p.Results.Warning, 'on')
        iswarning = 1;
else
        iswarning = 0;
end
switch p.Results.SendEL2
        case 'Text'
                option2do = 1;
        case 'EEG'
                option2do = 2;
        case 'EEG&Text'
                option2do = 3;
        case 'Workspace'
                option2do = 4;
        case 'Workspace&Text'
                option2do = 5;
        case 'Workspace&EEG'
                option2do = 6;
        case {'All', 'Workspace&EEG&Text', 'Workspace&Text&EEG', 'EEG&Text&Workspace', 'EEG&Workspace&Text', 'Text&Workspace&EEG', 'Text&EEG&Workspace'}
                option2do = 7;
        otherwise
                %option2do = 0;
                error('ERPLAB says: You must specify what to do with the updated EVENTLIST.')
end
if isempty(file2) || strcmpi(file2,'no') || strcmpi(file2,'none')
        if isfield(EEG, 'EVENTLIST')
                if isfield(EEG.EVENTLIST, 'eventinfo')
                        if isempty(EEG.EVENTLIST(indexEL).eventinfo)
                                msgboxText = ['ERPLAB says: EVENTLIST.eventinfo structure is empty!\n'...
                                        'Use Create EVENTLIST before BINLISTER'];
                                error(sprintf(msgboxText));
                        end
                else
                        msgboxText =  ['ERPLAB says: EVENTLIST.eventinfo structure was not found!\n'...
                                'Use Create EVENTLIST before BINLISTER'];
                        error(sprintf(msgboxText));
                end
        else
                msgboxText =  ['ERPLAB says: EVENTLIST structure was not found!\n'...
                        'Use Create EVENTLIST before BINLISTER'];
                error(sprintf(msgboxText));
        end
        if ~isempty(EEG.EVENTLIST(indexEL)) && iswarning==1
                binaux    = [EEG.EVENTLIST(indexEL).eventinfo.bini];
                binhunter = binaux(binaux>0); %8/19/2009

                if ~isempty(binhunter) && ismember_bc2(option2do, [2 3 6 7])
                        msgboxText =  ['This dataset already has assigned bins.\n'...;
                                'Would you like to overwrite these bins?'];
                        title = 'ERPLAB: WARNING';
                        button =askquest(sprintf(msgboxText), title);
                        if strcmpi(button,'no')
                                disp('User canceled')
                                return
                        end
                end
        end

        %
        % Check for alphanumeric event codes
        %
        codelist = {EEG.EVENTLIST(indexEL).eventinfo.code};
        if nnz(cellfun(@ischar, codelist))>0
                msgboxText = ['Your dataset still has alphanumeric/string codes.\n\n'...
                        'You may try either eliminating nonnumeric information (see Create EVENTLIST)\n'...
                        'or remapping your events (see Create EVENTLIST Advanced)'];
                %title = 'ERPLAB: BDF Parsing Error';
                %errorfound(sprintf(msgboxText), title);
                %return
                error(msgboxText);
        end
        % save original EVENTLIST
        %ELaux = EEG.EVENTLIST(indexEL); % store original EVENTLIST

        %
        %  Reset Flags?
        %
        if flagrst==1
                % reset artifact flags
                EEG = resetflag(EEG, 255, indexEL);
        elseif flagrst==2
                % reset user flags
                EEG = resetflag(EEG, 65280, indexEL);
        elseif flagrst==3
                % reset ALL flags
                EEG = resetflag(EEG, indexEL);
        end
end

forbiddenCodeArray = p.Results.Forbidden;
ignoreCodeArray    = p.Results.Ignore;

if strcmpi(p.Results.UpdateEEG, 'on')
        updevent = 1;
else
        updevent = 0;
end
if strcmpi(p.Results.Report, 'on')
        reportable = 1;
else
        reportable = 0;
end
if ismember_bc2({p.Results.Saveas}, {'on','yes'})
        issaveas    = 1;
else
        issaveas    = 0;
end
if strcmpi(outputv, 'EVENTLIST') || strcmpi(outputv, 'EL')
        Voutput = 1; % EVENTLIST
else
        Voutput = 0; % EEG
end
if strcmpi(p.Results.History,'implicit')
        shist = 3; % implicit
elseif strcmpi(p.Results.History,'script')
        shist = 2; % script
elseif strcmpi(p.Results.History,'gui')
        shist = 1; % gui
else
        shist = 0; % off
end

%
% subroutine
%
%
% Call (neo)binlister
%
if reportable % workaround for a faster binlister... 02/28/12 JLC
        [EEG, EVENTLIST, binofbins, isparsenum] = binlister(EEG, file1, file2, file3, forbiddenCodeArray, ignoreCodeArray, reportable);
else
        [EEG, EVENTLIST, binofbins, isparsenum] = neobinlister2_fast(EEG, file1, file2, file3, forbiddenCodeArray, ignoreCodeArray, 0);
end
if isparsenum==0
        % parsing was not approved
        msgboxText = ['Bin descriptor file contains errors!\n'...
                'For details, please read command window messages.'];
        title = 'ERPLAB: BDF Parsing Error';
        errorfound(sprintf(msgboxText), title);
        %EEG.EVENTLIST(indexEL) = ELaux;
        return
end
if nnz(binofbins)>=1
        if ~isempty(EVENTLIST)
                if Voutput==1 % EVENTLIST is the output only (for scripting)
                        clear EEG
                        EEG = EVENTLIST;
                        return
                end
                if iseegstruct(EEG) && ismember_bc2(option2do, [2 3 6 7])
                        EEG =  pasteeventlist(EEG, EVENTLIST, 1, indexEL);
                        if updevent && issaveas
                                EEG = pop_overwritevent(EEG, 'code', 'History', 'off');
                        elseif updevent && ~issaveas
                                EEG = pop_overwritevent(EEG, 'code', 'History', 'off');
                        end
                end
                if ismember_bc2(option2do, [4 5 6 7])
                        assignin('base','EVENTLIST',EVENTLIST);  % send EVENTLIST structure to WORKSPACE, August 22, 2008
                        disp('EVENTLIST structure was sent to WORKSPACE.')
                end
                if ismember_bc2(option2do, [1 3 5 7]) && ~isempty(file2) && ~strcmpi(file2,'no') && ~strcmpi(file2,'none')
                        disp('A text file version of your EVENTLIST was created.')
                end

                %
                % Generate equivalent command (for history)
                %
                skipfields = {'EEG', 'ERP', 'Saveas', 'Warning', 'History'};
                fn  = fieldnames(p.Results);
                com = sprintf( '%s  = %s( %s ', inputname(1), mfilename, inputname(1));
                for q=1:length(fn)
                        fn2com = fn{q};
                        if ~ismember_bc2(fn2com, skipfields)
                                fn2res = p.Results.(fn2com);
                                if ~isempty(fn2res)
                                        if ischar(fn2res)
                                                if ~strcmpi(fn2res,'off') && ~strcmpi(fn2res,'no') && ~strcmpi(fn2res,'none')
                                                        com = sprintf( '%s, ''%s'', ''%s''', com, fn2com, fn2res);
                                                end
                                        else
                                                if iscell(fn2res)
                                                        fn2resstr = vect2colon(cell2mat(fn2res), 'Sort','on');
                                                        fnformat = '{%s}';
                                                else
                                                        fn2resstr = vect2colon(fn2res, 'Sort','on');
                                                        fnformat = '%s';
                                                end
                                                com = sprintf( ['%s, ''%s'', ' fnformat], com, fn2com, fn2resstr);
                                        end
                                end
                        end
                end
                com = sprintf( '%s );', com);
        else
                msgboxText = ['Bin descriptor file contains errors!\n'...
                        'For details, please read command window messages.'];
                title = 'ERPLAB: BDF Parsing Error';
                errorfound(sprintf(msgboxText), title);
                EEG.EVENTLIST(indexEL) = ELaux;
                return
        end
else
        msgboxText =  ['Bins were not found!\n'...
                'Try with other BDF or modify the current one.'];
        title = 'ERPLAB: Binlister Error';
        errorfound(sprintf(msgboxText), title);
        disp('binlister process was cancel.')
        %EEG.EVENTLIST(indexEL) = ELaux;
        return
end
% get history from script. EEG
switch shist
        case 1 % from GUI
                com = sprintf('%s %% GUI: %s', com, datestr(now));
                %fprintf('%%Equivalent command:\n%s\n\n', com);
                displayEquiComERP(com);
        case 2 % from script
                EEG = erphistory(EEG, [], com, 1);
        case 3
                % implicit
        otherwise %off or none
                com = '';
end

%
% Completion statement
%
msg2end

if issaveas && ismember_bc2(option2do, [2 3 6 7]) && iseegstruct(EEG)
        [ALLEEG, EEG, CURRENTSET] = pop_newset( ALLEEG, EEG, CURRENTSET);
end
if exist('ALLEEG', 'var')
        if ~isfield(ALLEEG,'data')
                [ALLEEG(1:length(ALLEEG)).data] = deal([]);
        end
end
return

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
% FitDataNDFast : Fits several ND Gaussians to N-dimensional image data
% Synopsis: [fitted,params,myfunct] = FitDataNDFast(InitParm,Data,numdims,maxiter,method)
% AFunctString should use the variable 'x' as a vector and x{1}, x{2} to access its components,
% and the values to fit: c(1),c(2),...
%
% InitParm : Vector of initial parameters first row is global parameters,
% meaning is as follows: [global background, global width, intensity, posy, posx, intensity2, posy2, posx2...]
% other rows are local parameters
% Data : Experimental multidimensional (e.g. image) data to be fitted
% maxiter : maximum number of iterations in fit (default = 300)
% method : figure of merit to use for fitting 'mes', 'idiv' or 'fidiv'
%
% Example:
% a=noise(51.2*exp(-((xx(20,20)-2).^2+(yy(20,20)-1.2).^2)/20)+33,'poisson')
% [params,res,fitted]=FitDataNDFast([30 15 40 2 2],a,2,300,'idiv')
% fitted is the residual image (depending on method)
% params contains the result of the fit

function [params,res,fitted,residual] = FitDataNDFast(InitParm,Data,numdims,maxiter,mymethod)

if nargin < 3
    numdims = 2;
end
if nargin < 4
    maxiter = 3000;
end

if nargin < 5
    mymethod = 'mse';
end

%if nargin < 5
%    placement = 'right';
%end

%fixedparams=[0 20];
MultiGaussMSE(double(Data),mymethod,numdims)

[res, fitted] = MultiGaussMSE(InitParm');
%params=fixedparams;

if (0)  % old method
options=optimset('Display','notify','TolX',10^-4,'TolFun',10^-8,'MaxIter',maxiter);
[params,msevalue]=fminsearch(@MultiGaussMSE,InitParm',options);
else    % new method with smarter amd faster optimisation routine from http://www.cs.ubc.ca/~schmidtm
options=struct('Display','off','notify',1,'numDiff',1,'TolX',10^-4,'TolFun',10^-8,'MaxIter',maxiter);
[params,msevalue,moreinfo]=minFunc(@MultiGaussMSE,InitParm',options);
end

[res, fitted,residual] = MultiGaussMSE(params);
params = params';
fitted=dip_image(fitted);
residual=dip_image(residual);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB