text:
function [ResultsOrganEss] = organEssentiality(model, LPSolver)
% This function computes the organ essentiality in a whole-body model by
% setting each organ's value to zero in the whole-body objective reaction,
% setting all organ-specific reaction bounds to zero (lower and upper),
% and then computes whether a non-zero flux through this objective is still
% possible.
%
% function [ResultsOrganEss] = organEssentiality(model, LPSolver)
%
% INPUT
% model             model structure (whole-body metabolic model)
% LPSolver          Define LPSolver ('ILOGcomplex';
%                   'tomlab_cplex' -default)
%
% OUTPUT
% ResultsOrganEss   Contains the maximally possible flux value for the
%                   whole-body reaction for each organ. Col 1: Organ name,
%                   Col 2: Max flux value, Col 3: Min flux value (if
%                   minimization is activated; default: inactive), Col 4:
%                   Solver status: 1 = feasible, 5 = feasible with
%                   numerical difficulties (rescaling issues), 3 =
%                   infeasible)
%
% Ines Thiele 2016
% Ines Thiele, added option to specify LPSolver - 10/2018

if ~exist('LPSolver','var')
    LPSolver = 'tomlab_cplex';
end

global useSolveCobraLPCPLEX
useSolveCobraLPCPLEX

[solverOK, solverInstalled] = changeCobraSolver(LPSolver, 'LP');

% reset the bounds on the whole-body objective
model = changeRxnBounds(model,'Whole_body_objective_rxn',0,'l');
model = changeRxnBounds(model,'Whole_body_objective_rxn',1000,'u');

% define whole-body objective as objective function
model = changeObjective(model,'Whole_body_objective_rxn');

sex = model.sex;

% List of organs
getOrganWeightFraction;
OrganLists;
for i  = 1 : length(OrgansListShort)
    modelOrganEss = model;
    % redefine objective by removing the organ from the whole-body reaction
    O = strmatch(OrgansListShort{i},ObjectiveComponents);
    R = find(ismember(modelOrganEss.rxns,'Whole_body_objective_rxn'));
    M = find(ismember(modelOrganEss.mets,strcat(ObjectiveComponents{O},'_dummy_objective')));
    if useSolveCobraLPCPLEX
        if ~isfield(modelOrganEss,'A')
            error('model.A missing')
        end
        modelOrganEss.A(M,R)=0; % no requirement of this objective part in OF
    else
        modelOrganEss.S(M,R)=0; % no requirement of this objective part in OF
    end
    modelOrganEss = changeObjective(modelOrganEss,'Whole_body_objective_rxn');
    % set all reaction bounds in this organ to 0
    R1 = strmatch(OrgansListShort{i},modelOrganEss.rxns);
    modelOrganEss.lb(R1)=0;
    modelOrganEss.ub(R1)=0;
    % maximize the whole-body reaction
    modelOrganEss.osenseStr = 'max';
    tic;
    if useSolveCobraLPCPLEX
        [FBA,~]=solveCobraLPCPLEX(modelOrganEss,1,0,0,[],0,LPSolver);
        FBA.f=FBA.obj;
        FBA.v=FBA.full;
    else
        FBA = optimizeWBModel(modelOrganEss);
    end
    timeTaken = toc;
    fprintf('%u%s%s%s%f\n',timeTaken,' sec. ',OrgansListShort{i},' obj = ',FBA.f)
    ResultsOrganEss(i,1)=OrgansListShort(i);

    if useSolveCobraLPCPLEX
        feasible = FBA.origStat == 1 || FBA.origStat == 5;
    else
        feasible = FBA.stat == 1;
    end
    if feasible
        ResultsOrganEss{i,2}=num2str(FBA.v(modelOrganEss.c~=0)); % max
        if 0 % also compute the minimal possible flux through the objective
            modelOrganEss.osenseStr = 'min';
            tic;
            if useSolveCobraLPCPLEX
                [FBA,LPProblem]=solveCobraLPCPLEX(modelOrganEss,1,0,0,[],0,LPSolver);
            else
                FBA = optimizeWBModel(modelOrganEss);
            end
            timeTaken = toc;
            ResultsOrganEss{i,3}=num2str(FBA.v(modelOrganEss.c)~=0);%min
        end
        %feasible
        ResultsOrganEss{i,4}=num2str(FBA.origStat);
    else
        if FBA.stat == -1 ||  FBA.stat == 2
            %display solution in case there is a problem
            fprinf('%s\n',['Problem with FBA for organ: ' ResultsOrganEss{i,1} '. FBA solution is:'])
            FBA
        end

        %infeasible or otherwise
        ResultsOrganEss{i,2}=NaN; % min
        ResultsOrganEss{i,3}=NaN;%max
        ResultsOrganEss{i,4}=num2str(FBA.origStat);
    end
end

author: opencobra

repo: cobratoolbox


text:
% CROSSF - Returns estimates and plot of event-related coherence (ERC) changes
%            between data from two input channels. The lower panel gives the
%            coherent phase difference between the processes. In this panel, for Ex.
%               -90 degrees (blue) means xdata leads ydata by a quarter cycle.
%                90 degrees (orange) means ydata leads xdata by a quarter cycle.
%            Click on each subplot to view separately and zoom in/out.
%
% Function description:
%            Uses EITHER fixed-window, zero-padded FFTs (faster) OR constant-Q
%            0-padded DFTs (better sensitivity), both Hanning-tapered. Output
%            frequency spacing is the lowest frequency (srate/winsize) divided
%            by the padratio.
%
%            If number of output arguments > 4, then bootstrap statistics are
%            computed (from a distribution of 200 (NACCU) surrogate baseline
%            data epochs) for the baseline epoch, and non-significant features
%            of the output plots are zeroed (e.g., plotted in green). Baseline
%            epoch is all windows with center times < 0 (MAX_BASELN)
%
%            If number of output arguments > 5, coherency angles (lags) at
%            significant coherency (time,frequency) points are plotted as well.
%
% Usage:
%      >> [coh,mcoh,timesout,freqsout,cohboot,cohangles] = crossf(xdata,ydata,...
%                                              frames,tlimits,titl,          ...
%                                              srate,cycles,winsize,timesout,...
%                                              padratio,maxfreq,alpha,verts);
%
% Inputs:
%       xdata       = first single-channel (1,frames*nepochs) data  {none}
%       ydata       = second single-channel (1,frames*nepochs) data {none}
%       frames      = frames per epoch                        {768}
%       tlimits     = epoch time limits (ms) [mintime maxtime]{-1000 2000}
%       titl        = figure title                            {none}
%       srate       = data sampling rate (Hz)                 {256}
%       cycles      = >0 -> number of cycles in each analysis window (slower)
%                     =0 -> use FFT (constant window length)  {0}
%       winsize     = cycles==0: data subwindow length (2^k<frames)
%                     cycles >0: *longest* window length to use;
%                     determines the lowest output frequency  {~frames/8}
%       timesout    = number of output times (int<frames-winsize){200}
%       padratio    = FFT-length/winsize (2^k)                {2}
%                     Multiplies the number of output frequencies.
%       maxfreq     = maximum frequency to plot (Hz)          {50}
%       alpha       = Two-tailed bootstrap signif. probability {0.02}
%                     Sets n.s. plotted output values to green (0).
%                     NOTE that it requires at least FIVE output arguments!
%       verts       = times of vertical lines (other than time 0) {none}
%       caxma       = color axis maximum (magnitude) {default: from data}
%
% Outputs:
%       coh         = between-channel coherency changes (nfreqs,timesout)
%       mcoh        = vector of mean baseline coherence at each frequency
%       timesout    = vector of output times (subwindow centers) in ms.
%       freqsout    = vector of frequency bin centers in Hz.
%       cohboot     = [2,nfreqs] matrix of [lower;upper] coh significance diffs.
%       cohangle    = coherency angles (nfreqs,timesout)
%
% Note: when cycles==0, nfreqs is total number of FFT frequencies.
%
% Authors: Sigurd Enghoff & Scott Makeig, SCCN/INC/UCSD, La Jolla, 1998
%
% See also: TIMEF

% Copyright (C) 8/1/98 Sigurd Enghoff & Scott Makeig, SCCN/INC/UCSD
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 11-20-98 defined LINEWIDTH constant -sm
% 04-01-99 made number of frequencies consistent -se
% 06-29-99 fixed constant-Q freq indexing -se
% 08-13-99 added cohangle plotting -sm
% 08-20-99 made bootstrap more efficient -sm
% 08-24-99 allow nan values introduced by possible EVENTLOCK preproc. -sm
% 03-16-00 added lead/lag interpretation to help msg - sm & eric visser
% 03-16-00 added AXCOPY feature -sm & tpj
% 04-20-00 fixed Rangle sign for wavelets, added verts array -sm
% 01-22-01 corrected help msg when nargin<2 -sm & arno delorme
% 01-25-02 reformated help & license, added links -ad

function [R,mbase,times,freqs,Rboot,Rangle,Rsignif] = crossf(X,Y,epoch,timelim,ftitle,Fs,varwin,winsize,nwin,oversmp,maxfreq,alpha,verts,caxmax)

% Constants set here:
MAX_BASELN      = 0;            % Windows with center times < this are in baseline.
NACCU           = 200;			% Number of sub-windows to accumulate
if nargin>13
  COH_CAXIS_LIMIT = caxmax;
else
  COH_CAXIS_LIMIT = 0;          % 0 -> use data limits; else positive value
end
                                % giving symmetric +/- caxis limits.
AXES_FONT       = 10;
LINEWIDTH       = 2;
TITLE_FONT      = 8;
ANGLEUNITS      = 'deg';        % angle plotting units - 'ms' or 'deg'

% Commandline arg defaults:
DEFAULT_EPOCH	= 768;			% Frames per epoch
DEFAULT_TIMELIM = [-1000 2000];	% Time range of epochs (ms)
DEFAULT_FS		= 256;			% Sampling frequency (Hz)
DEFAULT_NWIN	= 200;			% Number of windows = horizontal resolution
DEFAULT_VARWIN	= 0;			% Fixed window length or base on cycles.
								% =0: fix window length to nwin
								% >0: set window length equal varwin cycles
								%     bounded above by winsize, also determines
								%     the min. freq. to be computed.
DEFAULT_OVERSMP	= 2;			% Number of times to oversample = vertical resolution
DEFAULT_MAXFREQ = 50;			% Maximum frequency to display (Hz)
DEFAULT_TITLE	= '';			% Figure title
DEFAULT_ALPHA   = 0.02;			% Default two-sided significance probability threshold
MARGIN          = 0.12;         % width of marginal plots
DEFAULT_VERTS   = [];           % default no vertical lines

if (nargin < 2)
	help crossf
	return
end

if (min(size(X))~=1 || length(X)<2)
	fprintf('crossf(): xdata must be a row or column vector.\n');
    return
elseif (min(size(Y))~=1 || length(Y)<2)
	fprintf('crossf(): ydata must be a row or column vector.\n');
    return
elseif (length(X) ~= length(Y))
	fprintf('crossf(): xdata and ydata must have same length.\n');
    return
end

if (nargin < 3)
	epoch = DEFAULT_EPOCH;
elseif (~isnumeric(epoch) || length(epoch)~=1 || epoch~=round(epoch))
	fprintf('crossf(): Value of frames must be an integer.\n');
    return
elseif (epoch <= 0)
	fprintf('crossf(): Value of frames must be positive.\n');
    return
elseif (rem(length(X),epoch) ~= 0)
	fprintf('crossf(): Length of data vectors must be divisible by frames.\n');
    return
end

if (nargin < 4)
	timelim = DEFAULT_TIMELIM;
elseif (~isnumeric(timelim) || sum(size(timelim))~=3)
	error('crossf(): Value of tlimits must be a vector containing two numbers.');
elseif (timelim(1) >= timelim(2))
	error('crossf(): tlimits interval must be [min,max].');
end

if (nargin < 5)
	ftitle = DEFAULT_TITLE;
elseif (~ischar(ftitle))
	error('crossf(): Plot title argument must be a quoted string.');
end

if (nargin < 6)
	Fs = DEFAULT_FS;
elseif (~isnumeric(Fs) || length(Fs)~=1)
	error('crossf(): Value of srate must be a number.');
elseif (Fs <= 0)
	error('crossf(): Value of srate must be positive.');
end

if (nargin < 7)
	varwin = DEFAULT_VARWIN;
elseif (~isnumeric(varwin) || length(varwin)~=1)
	error('crossf(): Value of cycles must be a number.');
elseif (varwin < MAX_BASELN)
	error('crossf(): Value of cycles must be either zero or positive.');
end

if (nargin < 8)
	winsize = max(pow2(nextpow2(epoch)-3),4);
elseif (~isnumeric(winsize) || length(winsize)~=1 || winsize~=round(winsize))
	error('crossf(): Value of winsize must be an integer number.');
elseif (winsize <= 0)
	error('crossf(): Value of winsize must be positive.');
elseif (varwin == 0 && pow2(nextpow2(winsize)) ~= winsize)
	error('crossf(): Value of winsize must be an integer power of two [1,2,4,8,16,...]');
elseif (winsize > epoch)
	error('crossf(): Value of winsize must be less than epoch length.');
end

if (nargin < 9)
	nwin = DEFAULT_NWIN;
elseif (~isnumeric(nwin) || length(nwin)~=1 || nwin~=round(nwin))
	error('crossf(): Value of nwin must be an integer number.');
elseif (nwin <= 0)
	error('crossf(): Value of nwin must be positive.');
end
if (nwin > epoch-winsize)
	error('crossf(): Value of nwin must be <= epoch-winsize.');
end

if (nargin < 10)
	oversmp = DEFAULT_OVERSMP;
elseif (~isnumeric(oversmp) || length(oversmp)~=1 || oversmp~=round(oversmp))
	error('crossf(): Value of oversmp must be an integer number.');
elseif (oversmp <= 0)
	error('crossf(): Value of oversmp must be positive.');
elseif (pow2(nextpow2(oversmp)) ~= oversmp)
	error('crossf(): Value of oversmp must be an integer power of two [1,2,4,8,16,...]');
end

if (nargin < 11)
	maxfreq = DEFAULT_MAXFREQ;
elseif (~isnumeric(maxfreq) || length(maxfreq)~=1)
	error('crossf(): Value of maxfreq must be a number.');
elseif (maxfreq <= 0)
	error('crossf(): Value of maxfreq must be positive.');
end

if (nargin < 12)
	alpha = DEFAULT_ALPHA;
elseif (~isnumeric(alpha) || length(alpha)~=1)
	error('crossf(): Value of alpha must be a number.');
elseif (round(NACCU*alpha) < 2 || alpha > .5)
	fprintf('crossf(): Value of alpha must be in the range (~0,0.5]');
    return
else
    if round(NACCU*alpha)<1,
      alpha = 1/NACCU;
	  fprintf(...
        'Using alpha = %0.3f. To decrease, must raise NACCU in source code.\n',...
                       alpha);
    end
end
if (nargin < 13)
   verts = DEFAULT_VERTS;
end

if (varwin == 0) % FFT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	freqs = Fs/winsize*[1:2/oversmp:winsize]/2;
	win = hanning(winsize);

	R = zeros(oversmp*winsize/2,nwin);
	RR = zeros(oversmp*winsize/2,nwin);
	Rboot = zeros(oversmp*winsize/2,NACCU);

else % wavelet DFT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	freqs = Fs*varwin/winsize*[2:2/oversmp:winsize]/2;
	win = dftfilt(winsize,maxfreq/Fs,varwin,oversmp,1);

	R  = zeros(size(win,2),nwin);
	RR = zeros(size(win,2),nwin);
	Rboot = zeros(size(win,2),NACCU);
end

wintime = 500*winsize/Fs;
times = timelim(1)+wintime:(timelim(2)-timelim(1)-2*wintime)/(nwin-1):timelim(2)-wintime;

baseln = find(times < 0);
dispf = find(freqs <= maxfreq);
stp = (epoch-winsize)/(nwin-1);
trials = length(X)/epoch;

fprintf('\nComputing the Event-Related Cross-Coherence image\n');
fprintf(' based on %d trials of %d frames sampled at %g Hz.\n',...
                       trials,epoch,Fs);
fprintf('Trial timebase is %d ms before to %d ms after the stimulus\n',...
                       timelim(1),timelim(2));
fprintf('The frequency range displayed is %g-%g Hz.\n',min(dispf),maxfreq);
if varwin==0
  fprintf('The data window size is %d samples (%g ms).\n',winsize,2*wintime);
  fprintf('The FFT length is %d samples\n',winsize*oversmp);
else
  fprintf('The window size is %d cycles.\n',varwin);
  fprintf('The maximum window size is %d samples (%g ms).\n',winsize,2*wintime);
end
fprintf('The window is applied %d times\n',nwin);
fprintf(' with an average step size of %g samples (%g ms).\n',...
                       stp,1000*stp/Fs);
fprintf('Results are oversampled %d times.\n',oversmp);
if nargout>4
  fprintf('Bootstrap confidence limits will be computed based on alpha = %g\n',...
              alpha);
else
  fprintf('Bootstrap confidence limits will NOT be computed.\n');
end
if nargout>5
  fprintf(['Coherence angles will be imaged in ',ANGLEUNITS,' and saved.\n']);
end

fprintf('\nProcessing trial (of %d):',trials);
firstboot = 1;
Rn=zeros(1,nwin);
X = X(:)'; % make X and Y column vectors
Y = Y(:)';
for t=1:trials,
	if (rem(t,10) == 0)
		fprintf(' %d',t);
	end
    if rem(t,120) == 0
        fprintf('\n');
    end

	for j=1:nwin, % for each time window
		tmpX = X([1:winsize]+floor((j-1)*stp)+(t-1)*epoch);
		tmpY = Y([1:winsize]+floor((j-1)*stp)+(t-1)*epoch);

        if ~any(isnan(tmpX))
		  tmpX = tmpX - mean(tmpX);
		  tmpY = tmpY - mean(tmpY);

		  if varwin == 0 % use FFTs
			tmpX = win .* tmpX(:);
			tmpY = win .* tmpY(:);
			tmpX = fft(tmpX,oversmp*winsize);
			tmpY = fft(tmpY,oversmp*winsize);
			tmpX = tmpX(2:oversmp*winsize/2+1);
			tmpY = tmpY(2:oversmp*winsize/2+1);
		  else
			tmpX = win' * tmpX(:);
			tmpY = win' * tmpY(:);
		  end

          if nargout > 4
           if firstboot==1
             tmpsX = repmat(nan,length(tmpX),nwin);
             tmpsY = repmat(nan,length(tmpY),nwin);
             firstboot = 0;
           end
           tmpsX(:,j) = tmpX;
           tmpsY(:,j) = tmpY;
          end

		  RR(:,j) = tmpX.*conj(tmpY) ./ (abs(tmpX).*abs(tmpY)); % complex coher.
	      R(:,j) = R(:,j) + RR(:,j);
          Rn(j) = Rn(j)+1;
        end % ~any(isnan())
	end % time window

	if (nargout > 4) % get NACCU bootstrap estimates for each trial
        j=1;
		while j<=NACCU
           s = ceil(rand([1 2])*nwin); % random ints [1,nwin]
           tmpX = tmpsX(:,s(1));
           tmpY = tmpsY(:,s(2));
           if ~any(isnan(tmpX)) && ~any(isnan(tmpY))
		      RR = tmpX.*conj(tmpY) ./ (abs(tmpX).*abs(tmpY)); % complex coher.
		      Rboot(:,j) = Rboot(:,j) + RR;
              j = j+1;
           end
        end
	end

end % t = trial

fprintf('\nNow plotting...\n');

Rangle = angle(R);
if varwin ~= 0
   Rangle = -Rangle; % make lead/lag the same for FFT and wavelet analysis
end
R = abs(R) ./ (ones(size(R,1),1)*Rn);               % coherence magnitude
Rraw = R;						% raw coherence values
mbase = mean(R(:,baseln)');     % mean baseline coherence magnitude
% R = R - repmat(mbase',[1 nwin]);% remove baseline mean

if (nargout > 4) % bootstrap
	i = round(NACCU*alpha);
	Rboot = abs(Rboot) / trials; % normalize bootstrap magnitude to [0,1]
	Rboot = sort(Rboot');
	Rsignif = mean(Rboot(NACCU-i+1:NACCU,:)); % significance levels for Rraw
%	Rboot = [mean(Rboot(1:i,:))-mbase ; mean(Rboot(NACCU-i+1:NACCU,:))-mbase];
 	Rboot = [mean(Rboot(1:i,:)) ; mean(Rboot(NACCU-i+1:NACCU,:))];
end % NOTE: above, mean ?????

set(gcf,'DefaultAxesFontSize',AXES_FONT)
colormap(jet(256));

pos = get(gca,'position'); % plot relative to current axes
q = [pos(1) pos(2) 0 0];
s = [pos(3) pos(4) pos(3) pos(4)];
axis('off')

if nargout>5 % image coherence lag as well as coherence magnitude
  ybase = 0.5;
  subheight = 0.4;
  MARGIN = MARGIN*0.75;
else
  ybase = 0.0;
  subheight = 0.9;
end

%
% Image the coherence [% perturbations]
%
RR = R;
if (nargout > 4) % zero out (and 'green out') nonsignif. R values
	RR(find((RR > repmat(Rboot(1,:)',[1 nwin])) ...
              & (RR < repmat(Rboot(2,:)',[1 nwin])))) = 0;
end
if (nargout > 5) % zero out nonsignif. Rraw values
	Rraw(find(repmat(Rsignif',[1,size(Rraw,2)])>=Rraw))=0;
end

if COH_CAXIS_LIMIT == 0
    coh_caxis = max(max(R(dispf,:)))*[-1 1];
else
    coh_caxis = COH_CAXIS_LIMIT*[-1 1];
end

h(6) = axes('Units','Normalized',...
               'Position',[MARGIN ybase+MARGIN 0.9-MARGIN subheight].*s+q);

map=hsv(300); % install circular color map - green=0, yellow, orng, red, violet = max
              %                                         cyan, blue, violet = min
map = flipud([map(251:end,:);map(1:250,:)]);
map(151,:) = map(151,:)*0.9; % tone down the (0=) green!
colormap(map);

imagesc(times,freqs(dispf),RR(dispf,:),coh_caxis); % plot the coherence image

set(h(6),'Units','Normalized',...
               'Position',[MARGIN ybase+MARGIN 0.9-MARGIN subheight].*s+q);
hold on
plot([0 0],[0 freqs(max(dispf))],'--m','LineWidth',LINEWIDTH)
for i=1:length(verts)
  plot([verts(i) verts(i)],[0 freqs(max(dispf))],'--m','LineWidth',LINEWIDTH);
end
hold off
set(h(6),'YTickLabel',[],'YTick',[])
set(h(6),'XTickLabel',[],'XTick',[])
title('Event-Related Coherence')
h(8) = axes('Position',[.92 ybase+MARGIN .05 subheight].*s+q);
cbar(h(8),151:300,[0 coh_caxis(2)]); % use only positive colors (gyorv)
%                                      for coherences
%
% Plot delta-mean min and max coherence at each time point on bottom of image
%
h(10) = axes('Units','Normalized','Position',[MARGIN ybase 0.9-MARGIN MARGIN].*s+q);
Emax = max(R(dispf,:)); % mean coherence at each time point
Emin = min(R(dispf,:)); % mean coherence at each time point
plot(times,Emax,'b');
hold on
plot(times,Emin,'b');
plot([times(1) times(length(times))],[0 0],'LineWidth',0.7);
plot([0 0],[-500 500],'--m','LineWidth',LINEWIDTH);
for i=1:length(verts)
  plot([verts(i) verts(i)],[-500 500],'--m','LineWidth',LINEWIDTH);
end
axis([min(times) max(times) 0 max(Emax)*1.2])
tick = get(h(10),'YTick');
set(h(10),'YTick',[tick(1) ; tick(length(tick))])
set(h(10),'YAxisLocation','right')
midpos = get(h(10),'Position');
if nargout<6
 xlabel('Time (ms)')
end
ylabel('coh.')

%
% Plot mean baseline coherence at each freq on left side of image
%

h(11) = axes('Units','Normalized','Position',[0 ybase+MARGIN MARGIN subheight].*s+q);
E = mbase(dispf); % baseline mean coherence at each frequency
if (nargout > 4) % plot bootstrap significance limits (base mean +/-)
	plot(freqs(dispf),E,'m','LineWidth',LINEWIDTH); % plot mbase
    hold on
	% plot(freqs(dispf),Rboot(:,dispf)+[E;E],'g','LineWidth',LINEWIDTH);
	plot(freqs(dispf),Rboot([1 2],dispf),'g','LineWidth',LINEWIDTH);
	plot(freqs(dispf),Rsignif(dispf),'k:','LineWidth',LINEWIDTH);
	axis([freqs(1) freqs(max(dispf)) 0 max([E Rsignif])*1.2]);
else             % plot marginal mean coherence only
	plot(freqs(dispf),E,'LineWidth',LINEWIDTH);
	% axis([freqs(1) freqs(max(dispf)) min(E)-max(E)/3 max(E)+max(E)/3]);
	if ~isnan(max(E))
	   axis([freqs(1) freqs(max(dispf)) 0 max(E)*1.2]);
	end;
end

tick = get(h(11),'YTick');
set(h(11),'YTick',[tick(1) ; tick(length(tick))])
set(h(11),'View',[90 90])
xlabel('Freq. (Hz)')
ylabel('coh.')

if (length(ftitle) > 0) % plot title
	axes('Position',pos,'Visible','Off');
	h(12) = text(-.05,1.01,ftitle);
	set(h(12),'VerticalAlignment','bottom')
	set(h(12),'HorizontalAlignment','left')
	set(h(12),'FontSize',TITLE_FONT)
end
%
% Plot coherence time lags in bottom panel
%
if nargout>5
   h(13) = axes('Units','Normalized','Position',[MARGIN MARGIN 0.9-MARGIN subheight].*s+q);
   if strcmp(ANGLEUNITS,'ms')  % convert to ms
     Rangle = (Rangle/(2*pi)).*repmat(1000./freqs(dispf)',1,length(times));
     maxangle = max(max(abs(Rangle)));
   else
     Rangle = Rangle*180/pi; % convert to degrees
     maxangle = 180; % use full-cycle plotting
   end
   Rangle(find(Rraw==0)) = 0; % set angle at non-signif coher points to 0

   imagesc(times,freqs(dispf),Rangle(dispf,:),[-maxangle maxangle]); % plot the
   hold on                                             % coherence phase angles
   plot([0 0],[0 freqs(max(dispf))],'--m','LineWidth',LINEWIDTH); % zero-time line
   for i=1:length(verts)
     plot([verts(i) verts(i)],[0 freqs(max(dispf))],'--m','LineWidth',LINEWIDTH);
   end

   pos13 = get(h(13),'Position');
   set(h(13),'Position',[pos13(1) pos13(2) midpos(3) pos13(4)]);
   ylabel('Freq. (Hz)')
   xlabel('Time (ms)')
   h(14)=axes('Position',[.92 MARGIN .05 subheight].*s+q);
   cbar(h(14),0,[-maxangle maxangle]); % two-sided colorbar

   if (length(ftitle) > 0) % plot title
	axes('Position',pos,'Visible','Off');
	h(13) = text(-.05,1.01,ftitle);
	set(h(13),'VerticalAlignment','bottom')
	set(h(13),'HorizontalAlignment','left')
	set(h(13),'FontSize',TITLE_FONT)
   end
 end
axcopy(gcf);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [] = plotWavefield(f, c, fs, T, thetaArrivalAngles, phiArrivalAngles, amplitudes)

if isscalar(f)
    f = f*ones(1,numel(thetaArrivalAngles));
end

if ~exist('amplitudes','var')
    amplitudes = ones(1,numel(thetaArrivalAngles));
end

% Grid-points for calculating the field
[x, y] = meshgrid(-0.5:0.02:0.5,-0.5:0.02:0.5);
nPoints = size(x,1);
xPoints = x(:)';
yPoints = y(:)';
nSamples = T*fs;
t = 0:1/fs:nSamples/fs-1/fs;

signal_field = 0;
for k = 1:numel(thetaArrivalAngles)
    doa = squeeze(steeringVector(xPoints, yPoints, f(k), c, thetaArrivalAngles(k), phiArrivalAngles(k)));
    signal = 10^(amplitudes(k)/20)*doa*exp(1j*2*pi*f(k)*t);
    signal_field = signal_field + signal;
end

figure(11);clf
set(gcf,'color','w')
colormap('spring')
cmap = [1 1 1]*0.5;
%plot(0.5*cos(0:pi/50:2*pi),0.5*sin(0:pi/50:2*pi),'Color',[1 1 1],'linewidth',1.5)
hold on
set(gca,'color',[0 0 0],'xcolor',cmap,'ycolor',cmap,'zcolor',cmap)
set(gca,'XTickLabel',[],'YTickLabel',[],'ZTickLabel',[])
set(gca,'XMinorGrid','on','YMinorGrid','on','ZMinorGrid','on','MinorGridColor',[1 1 1],'MinorGridLineStyle','-')
axis([-0.5 0.5 -0.5 0.5 -0.5 0.5])
view(110,25)
text(0.46,-0.5,-0.45,'x','color',[1 1 1])
text(-0.5,0.46,-0.45,'y','color',[1 1 1])
text(-0.5,-0.5,0.46,'z','color',[1 1 1])

z = reshape(real(signal_field(:,1))'/25,nPoints,nPoints);
h = surf(x,y,z,'edgecolor','none','FaceAlpha',0.7');
h.ZDataSource = 'z';

title(['t = ' num2str(sprintf('%0.1f',t(1)*1e3)) ' ms'],'fontweight','normal')

for sample = 2:nSamples
    pause(0.01)

    z = reshape(real(signal_field(:,sample))'/25,nPoints,nPoints);
    refreshdata(h, 'caller')

    title(['t = ' num2str(sprintf('%0.1f',t(sample)*1e3)) ' ms'],'fontweight','normal')

end

author: jorgengrythe

repo: beamforming


text:
function pot = convert_to_pot(CPD, pot_type, domain, evidence)
% CONVERT_TO_POT Convert a softmax CPD to a potential
% pots = convert_to_pot(CPD, pot_type, domain, evidence)
%
% pots = CPD evaluated using evidence(domain)

ncases = size(domain,2);
assert(ncases==1); % not yet vectorized

sz = dom_sizes(CPD);
ns = zeros(1, max(domain));
ns(domain) = sz;

odom = domain(~isemptycell(evidence(domain)));
T = convert_to_table(CPD, domain, evidence);

switch pot_type
 case 'u',
  pot = upot(domain, sz, T, 0*myones(sz));
 case 'd',
  ns(odom) = 1;
  pot = dpot(domain, ns(domain), T);

 case {'c','g'},
  % Since we want the output to be a Gaussian, the whole family must be observed.
  % In other words, the potential is really just a constant.
  p = T;
  %p = prob_node(CPD, evidence(domain(end)), evidence(domain(1:end-1)));
  ns(domain) = 0;
  pot = cpot(domain, ns(domain), log(p));

 case 'cg',
  T = T(:);
  ns(odom) = 1;
  can = cell(1, length(T));
  for i=1:length(T)
    can{i} = cpot([], [], log(T(i)));
  end
  ps = domain(1:end-1);
  dps = ps(CPD.dpndx);
  cps = ps(CPD.cpndx);
  ddom = [dps CPD.self];
  cdom = cps;
  pot = cgpot(ddom, cdom, ns, can);

 case 'scg'
  T = T(:);
  ns(odom) = 1;
  pot_array = cell(1, length(T));
  for i=1:length(T)
    pot_array{i} = scgcpot([], [], T(i));
  end
  pot = scgpot(domain, [], [], ns, pot_array);

 otherwise,
  error(['unrecognized pot type ' pot_type])
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% SHIFT_TEST_1D.M
% Plot how step responses vary at level 1 to 4 with 16 shifted steps
% for DT CWT and DWT.
%
% Nick Kingsbury, Cambridge University, May 2002.

% Generate 16 step functions as columns of matrix X.
n = 120;
x = cumsum([zeros(n,16);diag(ones(16,1));zeros(n,16)]) - 0.5;

% Specify DT CWT filters here.
biort = 'near_sym_b';
qshift = 'qshift_b';

% Forward DT CWT on each column of x.
[Yl,Yh] = dtwavexfm(x,4,biort,qshift);

% Inverse DT CWT, one subband at a time, using gain_mask to select subbands.
z1 = dtwaveifm(Yl*0,Yh,biort,qshift,[1 0 0 0]);
z01 = dtwaveifm(Yl*0,Yh,biort,qshift,[0 1 0 0]);
z001 = dtwaveifm(Yl*0,Yh,biort,qshift,[0 0 1 0]);
z0001 = dtwaveifm(Yl*0,Yh,biort,qshift,[0 0 0 1]);
z0000 = dtwaveifm(Yl,Yh,biort,qshift,[0 0 0 0]);

% Check for perfect reconstruction: abs(error) < 1e-12.
z = z1 + z01 + z001 + z0001 + z0000;
DTCWT_error = max(abs(z(:)-x(:)))

% Plot 16 input steps, x, and 16 step responses at each level, z1 to z0000.
setfig(1)
set(gcf,'DefaultTextFontSize',12,'Color',[1 1 1]);
set(gcf,'numbertitle','off','name',['Step responses at 16 adjacent shifts']);
subplot('position',[0.15 0.1 0.3 0.8]);
sc = 1.2; % Scale offsets so sets of curves do not overlap.
on = ones(size(x,2)-1,1)/5;
offset = sc*cumsum([0;on;1;on;1;on;1;on;1;on;1;on]/4).'; % Offsets for each curve.
zdtcwt = [x+0.5 z1 z01 z001 z0001 z0000-0.5] - ones(size(x,1),1)*offset;
tt = -40:40;  % Limits of plot horizontally.
plot(tt,zdtcwt(n+8+tt,:),'-b');
axis off
text(0,-7*sc,'(a) Dual Tree CWT','horiz','c')
xpos = -42; % Position of text labels.
text(xpos,0*sc,'Input','horiz','r','vert','m');
text(xpos,-0.9*sc,'Wavelets','horiz','r','vert','m');
text(xpos,-1.4*sc,'Level 1','horiz','r','vert','m');
text(xpos,-2.4*sc,'Level 2','horiz','r','vert','m');
text(xpos,-3.4*sc,'Level 3','horiz','r','vert','m');
text(xpos,-4.4*sc,'Level 4','horiz','r','vert','m');
text(xpos,-5.5*sc,'Scaling fn','horiz','r','vert','m');
text(xpos,-6*sc,'Level 4','horiz','r','vert','m');
drawnow

% Now do DWT for comparison.

% Specify DWT filters here.
biort = 'antonini';

% Forward DWT on each column of x.
[Yl,Yh] = wavexfm(x,4,biort);

% Inverse DWT, one subband at a time, using gain_mask to select subbands.
z1 = waveifm(Yl*0,Yh,biort,[1 0 0 0]);
z01 = waveifm(Yl*0,Yh,biort,[0 1 0 0]);
z001 = waveifm(Yl*0,Yh,biort,[0 0 1 0]);
z0001 = waveifm(Yl*0,Yh,biort,[0 0 0 1]);
z0000 = waveifm(Yl,Yh,biort,[0 0 0 0]);

% Check for perfect reconstruction: abs(error) < 1e-12.
z = z1 + z01 + z001 + z0001 + z0000;
DWT_error = max(abs(z(:)-x(:)))

subplot('position',[0.47 0.1 0.3 0.8]);
on = ones(size(x,2)-1,1)/5;
offset = sc*cumsum([0;on;1;on;1;on;1;on;1;on;1;on]/4).';
zdwt = [x+0.5 z1 z01 z001 z0001 z0000-0.5] - ones(size(x,1),1)*offset;
tt = -40:40;
plot(tt,zdwt(n+8+tt,:),'-r');
axis off
% xlabel('(b) Real DWT')
text(0,-7*sc,'(b) Real DWT','horiz','c')

return

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
% Elmore delay sizing for a straight wire (GP)
% Boyd, "Problems in VLSI design" (Lecture)
% Written for CVX by Almir Mutapcic 02/08/06
%
% We consider the problem of finding optimal width profile
% for a straight wire segmented into N parts. We want to
% minimize the Elmore delay, subject to limits on wire width
% and the total area. We use a pi-model for each wire segment.
% Problem can be formulated as GP:
%
%   minimize   D
%       s.t.   w_min <= w <= w_max
%              area  <= Amax
%
% where variables are widths w (and arrival times T that are used
% to formulate the overall delay D expression).
%
% Important: We label root node as 1, and all the other nodes as
%            node_label_in_the_paper + 1 (due to Matlab's convention).
%            Also label nodes with increasing numbers downstream.

%********************************************************************
% user supplied data (problem constants and tree topology)
%********************************************************************
N = 10+1; % number of segments (including the root node which is labeled as 1)

% parent node array for the straight wire
% specifies which node is a unique parent for node i (always have a tree)
parent = [0:N-1];

% problem constants
Rsource = 0.1;
l = 1*ones(N-1,1);
alpha = 1*ones(N-1,1);
beta  = 1*ones(N-1,1);
gamma = 1*ones(N-1,1);

% load capacitance at each node
Cload = [0; ones(N-1,1)];

% minimum and maximum width and area specification
Wmin = 1;
Wmax = 10;
Amax = 50;

%********************************************************************
% derived data (computed from user's data)
%********************************************************************
% compute children cell array (evaluate who are children for each node)
children = cell(N,1);
leafs = [];
for node = [1:N]
  children{node} = find(parent == node);
  if isempty(children{node})
    leafs(end+1) = node; % leafs have no children
  end
end

%********************************************************************
% optimization
%********************************************************************

disp('Generating the tradeoff curve...')
Darray = []; widths = [];
for Amax = [10.05 10.5 11 12:2:20 22.5 25:5:60]
    fprintf( 'Amax = %5.2f: ', Amax );
    cvx_begin gp quiet
        % optimization variables
        variable w(N-1)     % wire width
        variable T(N)       % arrival time (Elmore delay to node i)

        % objective is the critical Elmore delay
        minimize( max( T(leafs) ) )
        subject to

          % wire segment resistance is inversely proportional to widths
          R = alpha.*l./w;
          R = [Rsource; R];

          % wire segment capacitance is an affine function of widths
          C_bar = beta.*l.*w + gamma.*l;
          C_bar = [0; C_bar];

          % compute common capacitances for each node (C_tilde in GP tutorial)
          C_tilde = cvx( zeros(N,1) );
          for node = [1:N]
            C_tilde(node,1) = Cload(node);
            for k = parent(node)
              if k > 0; C_tilde(node,1) = C_tilde(node,1) + C_bar(k); end;
            end
            for k = children{node}
              C_tilde(node,1) = C_tilde(node,1) + C_bar(k);
            end
          end

          % now compute total downstream capacitances
          C_total = C_tilde;
          for node = N:-1:1
            for k = children{node}
              C_total(node,1) = C_total(node,1) + C_total(k,1);
            end
          end

        % generate Elmore delay constraints
        R(1)*C_total(1) <= T(1,1);
        for node = 2:N
          R(node)*C_total(node) + T(parent(node),1) <= T(node,1);
        end

        % collect all the constraints
        sum(w.*l) <= Amax;
        Wmin <= w <= Wmax;
    cvx_end
    % display and store computed values
    fprintf('delay = %3.2f\n',cvx_optval);
    Darray = [Darray cvx_optval];
    widths = [widths w];
end

% indices of four taper designs on the tradeoff curve
Amax = [10.05 10.5 11 12:2:20 22.5 25:5:60];
A11ind = find(Amax == 11);
A20ind = find(Amax == 20);
A35ind = find(Amax == 35);
A50ind = find(Amax == 50);

% plot the tradeoff curve
figure, clf
plot(Darray,Amax, ...
     Darray(A11ind),Amax(A11ind),'ro',...
     Darray(A20ind),Amax(A20ind),'ro',...
     Darray(A35ind),Amax(A35ind),'ro',...
     Darray(A50ind),Amax(A50ind),'ro');
xlabel('Elmore delay D'); ylabel('Amax');
disp('Optimal tradeoff curve plotted.')

% plot four taper designs
figure, clf
w1 = widths(:,A50ind);
w2 = widths(:,A35ind);
w3 = widths(:,A20ind);
w4 = widths(:,A11ind);
plot_four_tapers(w1,w2,w3,w4);

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [imo,mask] = rbfwarp2d( im, ps, pd, varargin )
% Radial base function/Thin-plate spline 2D image warping.
% [imo,mask] = rbfwarp2d( im, ps, pd, method)
%   input:
%       im: image 2d matrix
%       ps: 2d source landmark [n*2]
%       pd: 2d destin landmark [n*2]
%       method:
%         'gau',r  - for Gaussian function   ko = exp(-|pi-pj|/r.^2);
%         'thin'   - for Thin plate function ko = (|pi-pj|^2) * log(|pi-pj|^2)
%   output:
%       imo  : output matrix
%       mask : mask for output matrix, 0/1 means out/in the border
%
%   Bookstein, F. L.
%   "Principal Warps: Thin Plate Splines and the Decomposition of Deformations."
%   IEEE Trans. Pattern Anal. Mach. Intell. 11, 567-585, 1989.
%
%   Code by WangLin
%   2015-11-5
%   wanglin193@hotmail.com

num_required_parameters = 3;
if nargin < num_required_parameters
    help rbfwarp2d.m
    return;
end

% initialize default parameters
[imh,imw,imc] = size(im);
r = 0.1*imw;
imo = zeros(imh,imw,imc);
% mask = zeros(imh,imw); % commented by Abdo

% parse parameters
if nargin > num_required_parameters
    iVarargin = 1;
    while iVarargin <= nargin - num_required_parameters
        switch lower(varargin{iVarargin})
            case 'thin'
                method = 't';
            case 'gau'
                method = 'g';
                r = varargin{iVarargin+1};
                iVarargin = iVarargin + 1;
        end
        iVarargin = iVarargin + 1;
    end
end

%% Training w with L
nump = size(pd,1);
num_center = size(ps,1);
K=zeros(nump,num_center);

for i=1:num_center
    %Inverse warping from destination!
    dx = ones(nump,1)*ps(i,:)-pd;
    K(:,i) = sum(dx.^2,2);
end

if( strcmpi(method,'g') )
    K = rbf(K,r);
elseif( strcmpi(method,'t') )
    K = ThinPlate(K);
end

% P = [1,xp,yp] where (xp,yp) are n landmark points (nx2)
P = [ones(num_center,1),pd];
% L = [ K  P;
%       P' 0 ]
L = [K,P;P',zeros(3,3)];
% Y = [x,y;
%      0,0]; (n+3)x2
Y = [ps;zeros(3,2)];
%w = inv(L)*Y;
w = L\Y;

%% Using w
[x,y] = meshgrid(1:imw,1:imh);
pt = [x(:), y(:)];

nump = size(pt,1);
Kp = zeros(nump,num_center);
for i=1:num_center
    dx = ones(nump,1)*ps(i,:)-pt;
    Kp(:,i) = sum(dx.^2,2);
end
if( strcmpi(method,'g') )
    Kp = rbf(Kp,r);
elseif( strcmpi(method,'t') )
    Kp = ThinPlate(Kp);
end

L = [Kp,ones(nump,1),pt];
ptall = L*w;

%reshape to 2d image
xd = reshape( ptall(:,1),imh,imw );
yd = reshape( ptall(:,2),imh,imw );

for i = 1:imc
    %tmpx = gpuArray(single(im(:,:,i)));
    imt= interp2( single(im(:,:,i)),xd,yd,'linear',0); %% 'linear' 'cubic' 'spline',    extrapval
%     imo(:,:,i) = uint8(imt);
    imo(:,:,i) = (imt);
end

mask = ~isnan(imt);
end

function ko = rbf(d,r)
    ko = exp(-d/r.^2);
end

function ko = ThinPlate(ri)
% k=(r^2) * log(r^2)
    r1i = ri;
    r1i((ri==0))=realmin; % Avoid log(0)=inf
    ko = (ri).*log(r1i);
end

author: AbdoKamel

repo: sidd-ground-truth-image-estimation


text:
function serialized_mat = serialize(mat)
% SYNTAX:
%   serialized_mat = serialize(mat)
% DESCRIPTION:
%   return the array of the serialized matrix

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

    serialized_mat = mat(:);
end

author: goGPS-Project


text:
% J. C. Spall, August 1999
% Written in support of text, Introduction to Stochastic Search and Optimization, 2003
% This program runs a GA with the tournament parent selection scheme
% used.  Elitism is included.  Parent selection includes the elitist
% elements.  The standard binary bit
% form is used here.  As usual, code works in terms of fitness
% values (higher better); results, however, are reported for the loss
% values of actual interest.  This code does not work with constraints
% on theta values other than those directly associated with thetamin,
% thetamax.
%
% This code calls the functions 'bit2num','num2bit',and'fitpop'.
%
% The user is expected to set a variable 'expect_fn' representing the
% expected number of function evaluations allowed.  The main "while"
% loop in the code tests for when the expected number of function
% evaluations exceeds 'expect_fn' and terminates the iteration on the
% first generation that would require a no. of function evaluations
% exceeding this value.  The test in the "while" loop uses a formula
% in ISSO, Chap. 10; the "identical" variable in the formula is replaced
% here by an actual running calculation of the identical total.  Note
% that the code actually makes more function evaluations than the
% expected total due to the desire to avoid excessive bookkeeping and
% slowing of the code.  However, the IDEA of terminating when the
% expecting running total exceeds a threshold is consistent with the
% idea of placing a premium on possibly expensive function evaluations.
%
% This code handles noisy loss functions via the fact that the variables
% 'loss' and 'lossfinal' can call different functions.  The global
% variable 'sigma' is used to generate the scaling for noise in the
% loss measurements.
%
clear all
global p N thetamin thetamax
cases=15;         %no. of Monte Carlo cases
N = 40;           % population size
p=10;             % dimension of theta
P_c =.6;          % crossover rate
P_m =0.0001;      % mutation rate
elite=2;                % no. of population elements for "elitism"
                        % (should pick s.t. N-elite = even no.)
expect_fn=4000;         % expected no. of function evaluations
M=4*ones(p,1);          % no. of places to right of decimal
thetamin=-1.6383*ones(p,1);  %lower bounds on theta elements
thetamax=1.6384*ones(p,1);   %upper bounds: ditto
thetamin_0=-1.6383*ones(p,1);   %lower bounds on theta for initialization
thetamax_0=1.6384*ones(p,1);    %upper bounds: ditto
lossfinaleval ='loss4thorder';    % choice of loss function for final perf.
                              % evaluation
loss='loss4thnoise';
lossfinalsq=0;          %variable for cum.(over 'cases')squared loss values
lossfinal=0;            %variable for cum. loss values
rand('seed',31415927);
randn('seed',3111113)
global z sigma;       %declaration of random var. used for normal noise
                      %generation in loss fn. calls given seed above;
                      %also sigma in noise (noise is dependent on theta)
sigma=0;            %multiplier in loss noise (lower bd. to s.d.)
fitness=zeros(N,1);                %dim.-setting statement
theta=zeros(p,1);                  %dim.-setting statement
% User warning if invalid N-elite value
if (round((N-elite)/2)~=(N-elite)/2)
   disp('WARNING: N-elite should be even number');
else
end
%
% Determination of no. bits/theta element
bits=zeros(p,1);
for i=1:p
   b=1;
   while 10^M(i)*(thetamax(i)-thetamin(i))>2^b-1
      b=b+1;
   end
   bits(i)=b;
end
length=ones(1,p)*bits;              %no. of bits/population element
thetabit=zeros(N,length);  %dummy statement setting dimensions of matrix
                           %containing bit values of all chromosomes
thetabit_noelite=zeros(N-elite,length);%dummy statement setting
                                       %dimension of matrix used to temp.
                                       %store offspring
parent1=zeros(1,length); %dummy dimension statement
parent2=zeros(1,length); %ditto
% Begin outer loop of 'cases' Monte Carlo runs
%
navg=0;
for c=1:cases
   c
   %
   % Initial Random Population
   %
   % Statements below allow intialization of the search using a random
   % placement in natural theta space.  Points are placed uniformly dist.
   % on the hypercube defined by thetamin_0, thetamax_0.  The points are
   % converted to a bit representation for use by the GA.
   %
   for i=1:N
      theta=(thetamax_0-thetamin_0).*rand(p,1)+thetamin_0;
      pointer=1;
      for j=1:p
         thetabit(i,pointer:pointer+bits(j)-1)=num2bit(theta(j),bits(j),...
            thetamin(j),thetamax(j));
         pointer=pointer+bits(j);
      end
   end

%thetabit = rand(N,length) > 0.5;  % sets array of 0/1 for
                                  % all N chromosomes
%loop below is included for testing purposes and artifical insertion of
%known solution (to ensure that alg. properly handles this)
%for i=1:bits(1):p*bits(1)
%   thetabit(N,i:i+bits(1)-1)=[0,1,1,1,1,1,1,1,1,1,1,1,1,1,1];
%end
%zer=bit2num(thetabit(N,1:1+bits(1)-1),bits(1),M(1),thetamin(1),thetamax(1))
%upper = zeros(n, 1);              % fills vector for reporting loss
%average = zeros(n, 1);            % ditto
%lower = zeros(n, 1);              % ditto
   fitness=fitpop(thetabit,bits,loss);%N by 1 vector of fitness values
   [bestfitness,bestindex]=max(fitness);  bestindex
   % Loop below produces the theta associated with the best fitness value
   pointer=1;
   for j=1:p
      theta(j)=bit2num(thetabit(bestindex,pointer:pointer+bits(j)-1),...
            bits(j),thetamin(j),thetamax(j));
      pointer=pointer+bits(j);
   end
   %
   % Main loop of GA
   %
   index=zeros(elite,1);
   %cumfit=zeros(N-elite,1);
   identical=0;   %initialization of variable that keeps track of no.
                  %of identical parents in selection pairs (due to
                  %resampling)
   n=1;
   while n< 1+(expect_fn-N+identical*(1-(1-P_m)^length))/(P_c*(N-elite)...
         +(1-P_c)*(N-elite)*(1-(1-P_m)^length))
      % Evaluate obj. function for each individual (an N-dim.vector)
      %fitness = fitpop(thetabit,bits,loss); %higher better
      % Invoke elitism for 'elite' best chromosomes
      temp_fit=fitness;
      for j=1:elite
         [junk,index(j)]=max(temp_fit);
         temp_fit(index(j))=min(temp_fit);
      end
      for j=1:2:N-elite-1
         % Tournament selection below.  Uses tournament size = 2.
         for k=1:2
            cand=ceil(N*rand(2,1));
            if fitness(cand(1))>fitness(cand(2))
               parent1_idx=cand(1);
            else
               parent1_idx=cand(2);
            end
            cand=ceil(N*rand(2,1));
            if fitness(cand(1))>fitness(cand(2))
               parent2_idx=cand(1);
            else
               parent2_idx=cand(2);
            end
         end
      % Crossover with above two parents
         if rand < P_c
            cross_point=ceil(rand*(length-1)); %an integer 1,...,length-1
            parent1=thetabit(parent1_idx,:);
            %above picks off the chrom. for parent1_idx
            parent2=thetabit(parent2_idx,:);
            identical=identical+2*xor(1,norm(parent1-parent2));
            %above picks off the chrom. for parent2_idx
            child1=[parent1(1:cross_point),...
               parent2(cross_point+1:length)];   %merges parents1&2
            child2=[parent2(1:cross_point),...
               parent1(cross_point+1:length)];   %merges parents1&2
            thetabit_noelite(j,:)=child1;%assigns child1 to new chrom. matrix
            thetabit_noelite(j+1,:)=child2;%ditto for child2
         else
            parent1=thetabit(parent1_idx,:);
               %above picks off the chrom. for parent1_idx
            parent2=thetabit(parent2_idx,:);
               %above picks off the chrom. for parent2_idx
            thetabit_noelite(j,:)=parent1;%assigns parent1 to chrom. matrix
            thetabit_noelite(j+1,:)=parent2;%ditto for parent2
         end
      end
      %
      % Mutation operator applied to thetabit_noelite vector (elite
      % chromosomes not subject to mutation)
      %
      mutate=rand(N-elite,length)<P_m;
      thetabit_noelite=xor(thetabit_noelite,mutate); %exclusive-or operator
                                               %invokes mutation
      % Put elite chromosomes back in population
      %
      for j=1:elite
        thetabit(j,:)=thetabit(index(j),:);
      end
      thetabit(elite+1:N,:)=thetabit_noelite;
      fitness=fitpop(thetabit,bits,loss);
      n=n+1;
   end
      [bestfitness,bestindex]=max(fitness);
   pointer=1;
   for j=1:p
      theta(j)=bit2num(thetabit(bestindex,pointer:pointer+bits(j)-1),...
      bits(j),thetamin(j),thetamax(j));
      pointer=pointer+bits(j);
   end
   lossvalue=feval(lossfinaleval,theta);lossvalue
   lossfinalsq=lossfinalsq+lossvalue^2;
   lossfinal=lossfinal+lossvalue;
   navg=navg+n-1;   	%counter for # of iters. taken
end
navg/cases
% standard dev. of loss values
if cases ~= 1
	disp('standard deviation of mean loss value')
   sd=((cases/(cases-1))^.5)*(lossfinalsq/cases-(lossfinal/cases)^2)^.5;
   sd=sd/(cases^.5)
else
end
% normalized loss values
disp('mean loss value over "cases" runs')
mean=lossfinal/cases

author: Sable

repo: mcbench-benchmarks


text:
function catalog = get_total(varargin)

%GET_TOTAL Read in the complete AEIC Total earthquake catalog
%  CATALOG = GET_TOTAL reads the entire "Total" database of AEIC
%  earthquakes. The output field CATALOG is a structure containing field
%  names that adhere fairly closely to their database definitions and
%  should be self explanatory. Origin time is given in Matlab date format.
%
% CATALOG = GET_TOTAL(DATABASE) specifies the path and name where
% the database can be found. If not specified the DATABASE name defaults to
% the location on the Seismology Linux Network. This function may work
% generically on other databases. However the table requirements and the
% preferred magnitude algorithm are tailoered specifically to the AEIC
% Total database.
%
% For each event GET_TOTAL distills the multiple possible
% magnitudes into a single creates a "preferred magnitude" based on the
% following order of preference:
%        1. Hrvd Mw
%        2. AEIC Mw
%        3. AEIC Ms
%        4. AEIC mb
%        5. AEIC ml
% Events with no magnitude are removed from the catalog. Md magnitudes are
% not currently considered (I believe all Md magnitudes also have an ml).

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

% CHECK ARGUMENTS
if ~admin.antelope_exists
    error('This function requires Antelope');
end
if numel(varargin)>0
    dbName = varargin{1};
else
    dbName = '/Seis/catalogs/aeic/Total/Total';
end;

% LOAD EVENTS FROM NETMAG TABLE WITH Mw MAGNITUDES
db = dbopen(dbName,'r');
db = dblookup(db,'','origin','','');
db1 = dblookup(db,'','netmag','','');
db = dbjoin(db,db1);
db = dbsubset(db,'magtype==''mw''');
[mw.orid,mw.magtype,mw.magnitude,mw.auth] = dbgetv(db,'orid','magtype','magnitude','netmag.auth');
dbclose(db);

% SELECT THE PREFERRED AUTH FOR Mw WHEN MULTIPLES EXIST
% Order hrvd before aeic
mw.auth = lower(mw.auth);
% order hrvd before aeic authored solutions
[tmp index] = sort(mw.auth);
index = flipud(index);
mw.orid = mw.orid(index);
mw.magtype = mw.magtype(index);
mw.magnitude = mw.magnitude(index);
mw.auth = mw.auth(index);
% remove duplicate Mw slecting prefered author
[tmp index tmp1] = unique(mw.orid);
mw.orid = mw.orid(index);
mw.magtype = mw.magtype(index);
mw.magnitude = mw.magnitude(index);
mw.auth = mw.auth(index);

% LOAD ORIGIN TABLE CATALOG
db = dbopen(dbName,'r');
db = dblookup(db,'','origin','','');
db = dbsort(db,'orid');
[catalog.time,catalog.lat,catalog.lon,catalog.depth,catalog.orid] = dbgetv(db,'time','lat','lon','depth','orid');
[catalog.nass,catalog.ndef,catalog.etype,catalog.auth] = dbgetv(db,'nass','ndef','etype','auth');
[catalog.ml,catalog.mb,catalog.ms] = dbgetv(db,'ml','mb','ms');
dbclose(db);
catalog.time = epoch2datenum(catalog.time);
f = find(catalog.ml==-999);     catalog.ml(f) = NaN;
f = find(catalog.mb==-999);     catalog.mb(f) = NaN;
f = find(catalog.ms==-999);     catalog.ms(f) = NaN;
catalog.mw = NaN * catalog.ml;

% ADD Mw TO ORIGIN CATALOG WHEN IT EXISTS
for n = 1:numel(mw.orid)
    f = find(catalog.orid==mw.orid(n));
    if numel(f)~=1
        error('orid mismatch');
    else
        catalog.mw(f) = mw.magnitude(n);
        %disp(['Updated orid ' num2str(catalog.orid(f)) '  ' num2str(mw.orid(n)) ]);
    end
end

% CREATE A PREFERRED MAGNITUDE
catalog.prefMagType = repmat({'none'},size(catalog.ml));
catalog.prefMagnitude = NaN * catalog.ml;
% set pref magnitude as ml
f = find(~isnan(catalog.ml));
catalog.prefMagType(f) = {'ml'};
catalog.prefMagnitude(f) = catalog.ml(f);
% set pref magnitude as mb
f = find(~isnan(catalog.mb));
catalog.prefMagType(f) = {'mb'};
catalog.prefMagnitude(f) = catalog.mb(f);
% set pref magnitude as ms
f = find(~isnan(catalog.ms));
catalog.prefMagType(f) = {'ms'};
catalog.prefMagnitude(f) = catalog.ms(f);
% set pref magnitude as mw
f = find(~isnan(catalog.mw));
catalog.prefMagType(f) = {'mw'};
catalog.prefMagnitude(f) = catalog.mw(f);

% REMOVE EVENTS WITH NO MAGNITUDE
f = find(~isnan(catalog.prefMagnitude));
numEventsAll = num2str(numel(catalog.prefMagnitude));
numEventsWithMag = num2str(numel(f));
catalog.time          = catalog.time(f);
catalog.lat           = catalog.lat(f);
catalog.lon           = catalog.lon(f);
catalog.depth         = catalog.depth(f);
catalog.orid          = catalog.orid(f);
catalog.nass          = catalog.nass(f);
catalog.ndef          = catalog.ndef(f);
catalog.etype         = catalog.etype(f);
catalog.auth          = catalog.auth(f);
catalog.ml            = catalog.ml(f);
catalog.mb            = catalog.mb(f);
catalog.ms            = catalog.ms(f);
catalog.mw            = catalog.mw(f);
catalog.prefMagType   = catalog.prefMagType(f);
catalog.prefMagnitude = catalog.prefMagnitude(f);

% STATUS LINE
dateEarliest = datestr(min(catalog.time),'yyyy/mm/dd HH:MM:SS');
dateMostRecent = datestr(max(catalog.time),'yyyy/mm/dd HH:MM:SS');
disp(['  Retrieved ' numEventsWithMag ' earthquakes with magnitudes (out of ' numEventsAll ' total)']);
disp(['       Earliest:    ' dateEarliest ]);
disp(['       Most recent: ' dateMostRecent ]);
disp('');

author: geoscience-community-codes

repo: GISMO


text:
function grid2mri_interp = grid_interp_mri_seeg(GridLoc, MRI)
% GRID_INTERP_MRI_SEEG: Interpolate a grid of points into a MRI.
%
% USAGE:  grid2mri_interp = grid_interp_mri(GridLoc, MRI)
%
% INPUT:
%     - GridLoc     : [Nx3] matrix, 3D positions of the volume grid points
%     - MRI         : Brainstorm MRI structure
% OUTPUT:
%     - grid2mri_interp : Sparse matrix [nVoxels, nGridLoc]

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2017-2020

% ===== CHECK MRI =====
% Check that MRI SCS is well defined
if ~isfield(MRI,'SCS') || ~isfield(MRI.SCS,'R') || ~isfield(MRI.SCS,'T') || isempty(MRI.SCS.R) || isempty(MRI.SCS.T)
    error(['MRI SCS (Subject Coordinate System) was not defined or subjectimage file is from another version of Brainstorm.' 10 10,...
           'Please define the SCS fiducials on this MRI.']);
end
cubeSize = size(MRI.Cube(:,:,:,1));
% Convert coordinates
GridLoc = cs_convert(MRI, 'scs', 'voxel', GridLoc);

% ===== CHECK VERTICES LOCATION =====
% Get all the vertices that are outside the MRI volume
iOutsideVert = find((GridLoc(:,1) >= cubeSize(1)) | (GridLoc(:,1) < 2) | ...
                    (GridLoc(:,2) >= cubeSize(2)) | (GridLoc(:,2) < 2) | ...
                    (GridLoc(:,3) >= cubeSize(3)) | (GridLoc(:,3) < 2));
% Compute percentage of vertices outside the MRI
percentOutside = length(iOutsideVert) / length(GridLoc);
% If more than 95% vertices are outside the MRI volume : exit with ar error message
if (percentOutside > .95)
    grid2mri_interp = [];
    java_dialog('error', ['SEEG/ECOG contacts are not registered with the MRI.' 10 'Please try to import the position of the contacts again.'], 'SEEG/ECOG -> MRI');
    return;
% If more than 40% vertices are outside the MRI volume : display warning message
elseif (percentOutside > .4)
    java_dialog('warning', 'SEEG/ECOG contacts do not seem to be registered with the MRI.', 'SEEG/ECOG -> MRI');
end

%% ===== COMPUTE INTERPOLATION =====
% Initialize interpolation matrix
N = size(GridLoc,1);
grid2mri_interp = sparse([],[],[], prod(cubeSize), N, N);
% Generate a 3x3 block
[Xs, Ys, Zs] = meshgrid([-1,0,1],[-1,0,1],[-1,0,1]);
% Add entry for each grid point
for i = 1:N
    % Find the closest voxel to the SEEG contact
    P = round(GridLoc(i,:));
    % Add the 3x3 block around the center
    X = min(max(1, P(1)+Xs(:)), cubeSize(1));
    Y = min(max(1, P(2)+Ys(:)), cubeSize(2));
    Z = min(max(1, P(3)+Zs(:)), cubeSize(3));
    indVol = sub2ind(cubeSize, X, Y, Z);
    % Set the voxel value to exactly the contact value
    grid2mri_interp(indVol,i) = 1;
end

author: brainstorm-tools


text:
% writecnt() - Write a Neuroscan continuous signal file.
%
% Usage:
%   >> writecnt(filename, CNT-dataset, varargin)
%
% Inputs:
%   filename - name of the file with extension
%   dataset -  name of the CNT-dataset, a structure with the following fields
%                   cntdataset.header
%                   cntdataset.electloc
%                   cntdataset.data
%                   cntdataset.Teeg
%                   cntdataset.event
%                   cntdataset.tag
%                   cntdataset.endtag
%
% Optional inputs:
%  't1'         - start at time t1, default 0
%  'sample1'    - start at sample1, default 0, overrides t1
%  'lddur'      - duration of segment to load, default = whole file
%  'ldnsamples' - number of samples to load, default = whole file,
%                 overrides lddur
%  'scale'      - ['on'|'off'] scale data to microvolt (default:'on')
%  'dataformat' - ['int16'|'int32'] default is 'int16' for 16-bit data.
%                 Use 'int32' for 32-bit data.
%
% Outputs:
%  file         - file containing NeuroScan CNT file with the continuous
%                 data and other information
%
% Authors:   Sean Fitzgibbon, Arnaud Delorme, Michiel Vestjens 2000-2004
%
% Known limitations:
%  For more see http://www.cnl.salk.edu/~arno/cntload/index.html

% Copyright (C) 2000 Sean Fitzgibbon, <psspf@id.psy.flinders.edu.au>
% Copyright (C) 2003 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function f = writecnt(filename,cntdataset,varargin)

if ~isempty(varargin)
	 WriteOptions=struct(varargin{:});
else WriteOptions = [];
end;

try, WriteOptions.t1;         catch, WriteOptions.t1=0; end
try, WriteOptions.sample1;    catch, WriteOptions.sample1=[]; end
try, WriteOptions.lddur;      catch, WriteOptions.lddur=[]; end
try, WriteOptions.ldnsamples; catch, WriteOptions.ldnsamples=[]; end
try, WriteOptions.scale;      catch, WriteOptions.scale='on'; end
try, WriteOptions.dataformat;  catch, WriteOptions.dataformat='int16'; end

sizeEvent1 = 8  ; %%% 8  bytes for Event1
sizeEvent2 = 19 ; %%% 19 bytes for Event2

type='cnt';
if nargin ==1
    scan=0;
end

h   = cntdataset.header;
e   = cntdataset.electloc;
dat = cntdataset.data;
eT  = cntdataset.Teeg;
ev2 = cntdataset.event;
t   = cntdataset.tag;
endtag=cntdataset.endtag;

fid = fopen(filename,'w');
disp(['Writing file ' filename ' ...'])

% HEADER : 900 bytes => Starts at 0h, finishes with 383h.

fwrite(fid,h.rev,'char');
fwrite(fid,h.nextfile,'long');
fwrite(fid,h.prevfile,'long');
fwrite(fid,h.type,'char');
fwrite(fid,h.id,'char');
fwrite(fid,h.oper,'char');
fwrite(fid,h.doctor,'char');
fwrite(fid,h.referral,'char');
fwrite(fid,h.hospital,'char');
fwrite(fid,h.patient,'char');
fwrite(fid,h.age,'short');
fwrite(fid,h.sex,'char');
fwrite(fid,h.hand,'char');
fwrite(fid,h.med,'char');
fwrite(fid,h.category,'char');
fwrite(fid,h.state,'char');
fwrite(fid,h.label,'char');
fwrite(fid,h.date,'char');
fwrite(fid,h.time,'char');
fwrite(fid,h.mean_age,'float');
fwrite(fid,h.stdev,'float');
fwrite(fid,h.n,'short');
fwrite(fid,h.compfile,'char');
fwrite(fid,h.spectwincomp,'float');
fwrite(fid,h.meanaccuracy,'float');
fwrite(fid,h.meanlatency,'float');
fwrite(fid,h.sortfile,'char');
fwrite(fid,h.numevents,'int');
fwrite(fid,h.compoper,'char');
fwrite(fid,h.avgmode,'char');
fwrite(fid,h.review,'char');
fwrite(fid,h.nsweeps,'ushort');
fwrite(fid,h.compsweeps,'ushort');
fwrite(fid,h.acceptcnt,'ushort');
fwrite(fid,h.rejectcnt,'ushort');
fwrite(fid,h.pnts,'ushort');
fwrite(fid,h.nchannels,'ushort');
fwrite(fid,h.avgupdate,'ushort');
fwrite(fid,h.domain,'char');
fwrite(fid,h.variance,'char');
fwrite(fid,h.rate,'ushort');
fwrite(fid,h.scale,'double');
fwrite(fid,h.veogcorrect,'char');
fwrite(fid,h.heogcorrect,'char');
fwrite(fid,h.aux1correct,'char');
fwrite(fid,h.aux2correct,'char');
fwrite(fid,h.veogtrig,'float');
fwrite(fid,h.heogtrig,'float');
fwrite(fid,h.aux1trig,'float');
fwrite(fid,h.aux2trig,'float');
fwrite(fid,h.heogchnl,'short');
fwrite(fid,h.veogchnl,'short');
fwrite(fid,h.aux1chnl,'short');
fwrite(fid,h.aux2chnl,'short');
fwrite(fid,h.veogdir,'char');
fwrite(fid,h.heogdir,'char');
fwrite(fid,h.aux1dir,'char');
fwrite(fid,h.aux2dir,'char');
fwrite(fid,h.veog_n,'short');
fwrite(fid,h.heog_n,'short');
fwrite(fid,h.aux1_n,'short');
fwrite(fid,h.aux2_n,'short');
fwrite(fid,h.veogmaxcnt,'short');
fwrite(fid,h.heogmaxcnt,'short');
fwrite(fid,h.aux1maxcnt,'short');
fwrite(fid,h.aux2maxcnt,'short');
fwrite(fid,h.veogmethod,'char');
fwrite(fid,h.heogmethod,'char');
fwrite(fid,h.aux1method,'char');
fwrite(fid,h.aux2method,'char');
fwrite(fid,h.ampsensitivity,'float');
fwrite(fid,h.lowpass,'char');
fwrite(fid,h.highpass,'char');
fwrite(fid,h.notch,'char');
fwrite(fid,h.autoclipadd,'char');
fwrite(fid,h.baseline,'char');
fwrite(fid,h.offstart,'float');
fwrite(fid,h.offstop,'float');
fwrite(fid,h.reject,'char');
fwrite(fid,h.rejstart,'float');
fwrite(fid,h.rejstop,'float');
fwrite(fid,h.rejmin,'float');
fwrite(fid,h.rejmax,'float');
fwrite(fid,h.trigtype,'char');
fwrite(fid,h.trigval,'float');
fwrite(fid,h.trigchnl,'char');
fwrite(fid,h.trigmask,'short');
fwrite(fid,h.trigisi,'float');
fwrite(fid,h.trigmin,'float');
fwrite(fid,h.trigmax,'float');
fwrite(fid,h.trigdir,'char');
fwrite(fid,h.autoscale,'char');
fwrite(fid,h.n2,'short');
fwrite(fid,h.dir,'char');
fwrite(fid,h.dispmin,'float');
fwrite(fid,h.dispmax,'float');
fwrite(fid,h.xmin,'float');
fwrite(fid,h.xmax,'float');
fwrite(fid,h.automin,'float');
fwrite(fid,h.automax,'float');
fwrite(fid,h.zmin,'float');
fwrite(fid,h.zmax,'float');
fwrite(fid,h.lowcut,'float');
fwrite(fid,h.highcut,'float');
fwrite(fid,h.common,'char');
fwrite(fid,h.savemode,'char');
fwrite(fid,h.manmode,'char');
fwrite(fid,h.ref,'char');
fwrite(fid,h.rectify,'char');
fwrite(fid,h.displayxmin,'float');
fwrite(fid,h.displayxmax,'float');
fwrite(fid,h.phase,'char');
fwrite(fid,h.screen,'char');
fwrite(fid,h.calmode,'short');
fwrite(fid,h.calmethod,'short');
fwrite(fid,h.calupdate,'short');
fwrite(fid,h.calbaseline,'short');
fwrite(fid,h.calsweeps,'short');
fwrite(fid,h.calattenuator,'float');
fwrite(fid,h.calpulsevolt,'float');
fwrite(fid,h.calpulsestart,'float');
fwrite(fid,h.calpulsestop,'float');
fwrite(fid,h.calfreq,'float');
fwrite(fid,h.taskfile,'char');
fwrite(fid,h.seqfile,'char');
fwrite(fid,h.spectmethod,'char');
fwrite(fid,h.spectscaling,'char');
fwrite(fid,h.spectwindow,'char');
fwrite(fid,h.spectwinlength,'float');
fwrite(fid,h.spectorder,'char');
fwrite(fid,h.notchfilter,'char');
fwrite(fid,h.headgain,'short');
fwrite(fid,h.additionalfiles,'int');
fwrite(fid,h.unused,'char');
fwrite(fid,h.fspstopmethod,'short');
fwrite(fid,h.fspstopmode,'short');
fwrite(fid,h.fspfvalue,'float');
fwrite(fid,h.fsppoint,'short');
fwrite(fid,h.fspblocksize,'short');
fwrite(fid,h.fspp1,'ushort');
fwrite(fid,h.fspp2,'ushort');
fwrite(fid,h.fspalpha,'float');
fwrite(fid,h.fspnoise,'float');
fwrite(fid,h.fspv1,'short');
fwrite(fid,h.montage,'char');
fwrite(fid,h.eventfile,'char');
fwrite(fid,h.fratio,'float');
fwrite(fid,h.minor_rev,'char');
fwrite(fid,h.eegupdate,'short');
fwrite(fid,h.compressed,'char');
fwrite(fid,h.xscale,'float');
fwrite(fid,h.yscale,'float');
fwrite(fid,h.xsize,'float');
fwrite(fid,h.ysize,'float');
fwrite(fid,h.acmode,'char');
fwrite(fid,h.commonchnl,'uchar');
fwrite(fid,h.xtics,'char');
fwrite(fid,h.xrange,'char');
fwrite(fid,h.ytics,'char');
fwrite(fid,h.yrange,'char');
fwrite(fid,h.xscalevalue,'float');
fwrite(fid,h.xscaleinterval,'float');
fwrite(fid,h.yscalevalue,'float');
fwrite(fid,h.yscaleinterval,'float');
fwrite(fid,h.scaletoolx1,'float');
fwrite(fid,h.scaletooly1,'float');
fwrite(fid,h.scaletoolx2,'float');
fwrite(fid,h.scaletooly2,'float');
fwrite(fid,h.port,'short');
fwrite(fid,h.numsamples,'ulong');
fwrite(fid,h.filterflag,'char');
fwrite(fid,h.lowcutoff,'float');
fwrite(fid,h.lowpoles,'short');
fwrite(fid,h.highcutoff,'float');
fwrite(fid,h.highpoles,'short');
fwrite(fid,h.filtertype,'char');
fwrite(fid,h.filterdomain,'char');
fwrite(fid,h.snrflag,'char');
fwrite(fid,h.coherenceflag,'char');
fwrite(fid,h.continuoustype,'char');
fwrite(fid,h.eventtablepos,'long');
fwrite(fid,h.continuousseconds,'float');
fwrite(fid,h.channeloffset,'long');
fwrite(fid,h.autocorrectflag,'char');
fwrite(fid,h.dcthreshold,'uchar');  % = 383H

% ELECT.DESCRIPTIONS : 75*n.channels bytes.
% Starts with 384h, finishes with (899+75*nchannels)dec
% 10 channels: 671h   % 12 channels: 707h
% 24 channels: A8Bh   % 32 channels: CE3h
% 64 channels: 1643h  % 128 channels: 2903h

for n = 1:h.nchannels
    fwrite(fid,e(n).lab,'char');
    fwrite(fid,e(n).reference,'char');
    fwrite(fid,e(n).skip,'char');
    fwrite(fid,e(n).reject,'char');
    fwrite(fid,e(n).display,'char');
    fwrite(fid,e(n).bad,'char');
    fwrite(fid,e(n).n,'ushort');
    fwrite(fid,e(n).avg_reference,'char');
    fwrite(fid,e(n).clipadd,'char');
    fwrite(fid,e(n).x_coord,'float');
    fwrite(fid,e(n).y_coord,'float');
    fwrite(fid,e(n).veog_wt,'float');
    fwrite(fid,e(n).veog_std,'float');
    fwrite(fid,e(n).snr,'float');
    fwrite(fid,e(n).heog_wt,'float');
    fwrite(fid,e(n).heog_std,'float');
    fwrite(fid,e(n).baseline,'short');
    fwrite(fid,e(n).filtered,'char');
    fwrite(fid,e(n).fsp,'char');
    fwrite(fid,e(n).aux1_wt,'float');
    fwrite(fid,e(n).aux1_std,'float');
    fwrite(fid,e(n).senstivity,'float');
    fwrite(fid,e(n).gain,'char');
    fwrite(fid,e(n).hipass,'char');
    fwrite(fid,e(n).lopass,'char');
    fwrite(fid,e(n).page,'uchar');
    fwrite(fid,e(n).size,'uchar');
    fwrite(fid,e(n).impedance,'uchar');
    fwrite(fid,e(n).physicalchnl,'uchar');
    fwrite(fid,e(n).rectify,'char');
    fwrite(fid,e(n).calib,'float');
end

% finding if 32-bits of 16-bits file
% ----------------------------------
begdata = ftell(fid);
enddata = h.eventtablepos;   % after data
if strcmpi(WriteOptions.dataformat, 'int16')
     nums    = (enddata-begdata)/h.nchannels/2;
else nums    = (enddata-begdata)/h.nchannels/4;
end;

% number of sample to write
% -------------------------
if ~isempty(WriteOptions.sample1)
   WriteOptions.t1      = WriteOptions.sample1/h.rate;
else
   WriteOptions.sample1 = WriteOptions.t1*h.rate;
end;
if strcmpi(WriteOptions.dataformat, 'int16')
     startpos = WriteOptions.t1*h.rate*2*h.nchannels;
else startpos = WriteOptions.t1*h.rate*4*h.nchannels;
end;
if isempty(WriteOptions.ldnsamples)
     if ~isempty(WriteOptions.lddur)
          WriteOptions.ldnsamples = round(WriteOptions.lddur*h.rate);
     else WriteOptions.ldnsamples = nums;
     end;
end;

% scaling data from microvolts
% ----------------------------
if strcmpi(WriteOptions.scale, 'on')
    disp('Scaling data .....')
    for i=1:h.nchannels
       bas=e(i).baseline;
       sen=e(i).senstivity;
       cal=e(i).calib;
       mf=sen*(cal/204.8);
       dat(i,:)=(dat(i,:)/mf)+bas;
   end
end

% write data
% ----------
disp('Writing data .....')
if type == 'cnt'
    channel_off = h.channeloffset/2;

      fseek(fid, startpos, 0);
	  if channel_off <= 1
          for temploop =1:WriteOptions.ldnsamples;
              fwrite(fid, dat(1:h.nchannels, temploop), WriteOptions.dataformat)';
          end;
 	  else
          for temploop =1:h.nchannels;
              fwrite(fid, dat(temploop, 1:channel_off), WriteOptions.dataformat)';
          end;

		  counter = 1;
 		  while counter*channel_off < WriteOptions.ldnsamples

              for temploop =1:h.nchannels;
                fwrite(fid, dat(temploop, counter*channel_off+1:counter*channel_off+channel_off), WriteOptions.dataformat)';
              end;
              counter = counter + 1;
		  end;
	  end;

      % write event table
      % -----------------
      frewind(fid);
      fseek(fid,h.eventtablepos,'bof');

      disp('Writing Event Table...')
      fwrite(fid,eT.teeg,'uchar');
      fwrite(fid,eT.size,'ulong');
      fwrite(fid,eT.offset,'ulong');

      if eT.teeg==2
          nevents=eT.size/sizeEvent2;
      elseif eT.teeg==1
          nevents=eT.size/sizeEvent1;
      else
          disp('No !!! teeg <> 2 and teeg <> 1');
          ev2 = [];
      end

     %%%% to change offset in points back to bytes
     if ~isempty(ev2)
          ev2p=ev2;
          ioff=900+(h.nchannels*75); %% initial offset : header + electordes desc
          for i=1:nevents
              ev2p(i).offset=((ev2p(i).offset + WriteOptions.sample1)*2*h.nchannels) +ioff; %% 2 short int end
          end
          ev2 = ev2p;
     end;

     if eT.teeg==2
          nevents=eT.size/sizeEvent2;
          for i=1:nevents
              ev2(i).stimtype      = fwrite(fid,ev2(i).stimtype,'ushort');
              ev2(i).keyboard      = fwrite(fid,ev2(i).keyboard,'char');
              ev2(i).keypad_accept = fwrite(fid,ev2(i).keypad_accept,'char');
              ev2(i).offset        = fwrite(fid,ev2(i).offset,'long');
              ev2(i).type          = fwrite(fid,ev2(i).type,'short');
              ev2(i).code          = fwrite(fid,ev2(i).code,'short');
              ev2(i).latency       = fwrite(fid,ev2(i).latency,'float');
              ev2(i).epochevent    = fwrite(fid,ev2(i).epochevent,'char');
              ev2(i).accept        = fwrite(fid,ev2(i).accept,'char');
              ev2(i).accuracy      = fwrite(fid,ev2(i).accuracy,'char');
          end
      elseif eT.teeg==1
          nevents=eT.size/sizeEvent1;
          for i=1:nevents
              ev2(i).stimtype      = fwrite(fid,ev2(i).stimtype,'ushort');
              ev2(i).keyboard      = fwrite(fid,ev2(i).keyboard,'char');
              ev2(i).keypad_accept = fwrite(fid,ev2(i).keypad_accept,'char');
              ev2(i).offset        = fwrite(fid,ev2(i).offset,'long');
          end;
      else
          disp('No !!! teeg <> 2 and teeg <> 1');
          ev2 = [];
      end
end

if size(endtag,1)>1
    fwrite(fid,endtag);
else
    fwrite(fid,t,'char');
end

fclose(fid);
disp(['Finished writing file ' filename ' ...'])

author: PatternRecognition

repo: OpenBMI


text:
function imsegs = im2superpixels(im, method, varargin )
    if nargin < 3
        % default parameters to generate superpixels
        switch method
            case 'pedro'
                sigma = 0.8;
                k = 100;
                min_size = 150;
            case 'SLIC'
                num_superpixel = 200;
            otherwise
                error( 'unknown method to generate superpixels.' );
        end
    else
        switch method
            case 'pedro'
                para = varargin{1};
                sigma = para(1);
                k = para(2);
                min_size = para(3);
            case 'SLIC'
                num_superpixel = varargin{1};
            otherwise
                error( 'unknown method to generate superpixels.' );
        end
    end

%     prefix = num2str(floor(rand(1)*10000000));
%     fn1 = ['./tmpim' prefix '.ppm'];
%     fn2 = ['./tmpimsp' prefix '.ppm'];
%     segcmd = ['E:\playerkk\code\MATLAB\segment\segment ', num2str(seg_para(1)),...
%         ' ', num2str(seg_para(2)), ' ', num2str(seg_para(3))];
%
%     imwrite(im, fn1);
%     system([segcmd ' ' fn1 ' ' fn2]);
    if isa(im, 'uint8')
        im = double(im);
    end

    if max(im(:)) < 10
        im = double(im * 255);
    end

    switch method
        case 'pedro'
            segim = mexSegment(im, sigma, k, int32(min_size));
        case 'SLIC'
            segim = uint8(mexSLIC(uint32(im), num_superpixel));
        otherwise
            error( 'unknown method to generate superpixels.' );
    end
    imsegs = processSuperpixelImage(segim);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function f = fgplvmSequenceObjective(xvec, model, Y, varargin)

% FGPLVMSEQUENCEOBJECTIVE Wrapper function for objective of a single sequence in latent space and the corresponding output sequence.
% FORMAT
% DESC provides a wrapper function for the negative log probability
% of a given data sequence under the posterior distribution of the
% Gaussian process induced by the training data..
% ARG vecx : time ordered locations in input space for the sequence
% placed as a vector using the matlab X(:)' notation.
% ARG model : the model structure for which the negative log
% probability of the given data under the posterior is to be computed.
% ARG Y : time ordered locations in data spaces for the sequence.
% ARG P1, P2, P3 ... : optional additional arguments to be passed
% to the model sequence log likelihood.
% RETURN f : the negative of the log probability of the given data
% sequence under the posterior distribution induced by the training data.
%
% SEEALSO : fgplvmCreate, fgplvmSequenceLogLikelihood, fgplvmOptimiseSequence
%
% COPYRIGHT : Neil D. Lawrence, 2006

% FGPLVM

X = reshape(xvec, size(Y, 1), model.q);
f = - fgplvmSequenceLogLikelihood(model, X, Y, varargin{:});

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function res = hlp_username
% Get the user's name.

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

res = char(java.lang.System.getProperty('user.name'));

repo: REST


text:
function HTotal=uHessian3D(xG,lRx,M)
%%UHESSIAN3D Determine the Hessian matrix (second derivative matrix) of a
%           a single direction cosine u component in 3D with respect to 3D
%           position. Relativity and atmospheric effects are not taken into
%           account. This function is essentially the same as uvHessian
%           without the v component.
%
%INPUTS: xG The 3XN target position vectors in the global coordinate system
%           with [x;y;z] components for which gradients are desired.
%       lRx The 3X1  position vector of the receiver. If omitted, the
%           receiver is placed at the origin.
%         M A 3X3 rotation matrix from the global Coordinate system to the
%           orientation of the coordinate system at the receiver. If
%           omitted, it is assumed to be the identity matrix.
%
%OUTPUTS: HTotal A 3X3XN matrix such that HTotal(:,:,i) is the Hessian
%          matrix with respect to the u component of the ith measurement.
%          The elements in the matrices for each component/ point are
%          ordered [d^2/(dxdx), d^2/(dxdy), d^2/(dxdz);
%                   d^2/(dydx), d^2/(dydy), d^2/(dydz);
%                   d^2/(dzdx), d^2/(dzdy), d^2/(dzdz)];
%          note that each matrix is symmetric (i.e.
%           d^2/(dydx)=d^2/(dxdy) ).
%
%A derivation of the components of the Jacobian is given in [1]. The
%Hessian is just one derivative higher.
%
%REFERENCES:
%[1] D. F. Crouse, "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems
%    Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%
%July 2021 David F.Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(M))
	M=eye(3,3);
end

if(nargin<2||isempty(lRx))
	lRx=zeros(3,1);
end

N=size(xG,2);
HTotal=zeros(3,3,1,N);
for curPoint=1:N
    %Convert the state into the local coordinate system.
    xLocal=M*(xG(1:3)-lRx(1:3));

    x=xLocal(1);
    y=xLocal(2);
    z=xLocal(3);

    r5=norm(xLocal)^5;

    H=zeros(3,3,1);

    %u Hessian values
    %du/(dxdx)
    H(1,1,1)=-((3*x*(y^2+z^2))/r5);
    %du/(dydy)
    H(2,2,1)=-((x*(x^2-2*y^2+z^2))/r5);
    %du/(dzdz)
    H(3,3,1)=-((x*(x^2+y^2-2*z^2))/r5);
    %du/(dxdy)
    H(1,2,1)=-((y*(-2*x^2+y^2+z^2))/r5);
    %du/(dydx)
    H(2,1,1)=H(1,2,1);
    %du/(dxdz)
    H(1,3,1)=-((z*(-2*x^2+y^2+z^2))/r5);
    %du/(dzdx)
    H(3,1,1)=H(1,3,1);
    %du/(dydz)
    H(2,3,1)=(3*x*y*z)/r5;
    %du/(dzdy)
    H(3,2,1)=H(2,3,1);

    %Rotate the values back into global coordinates
    HTotal(:,:,1,curPoint)=M'*H(:,:,1)*M;
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function b = triangulation_isequal(face,faced)

% triangulation_isequal - true if two triangulations are equals
%
%   b = triangulation_isequal(face,faced)
%
%   Copyright (c) 2008 Gabriel Peyre

I = setdiff(sort(face)',sort(faced)', 'rows');
b = isempty(I);

author: gpeyre

repo: matlab-toolboxes


text:
function FileString = generateFieldDescriptionFile(FileName)
% Generates the ModelFields.md file describing the Required and Optional Fields of a COBRA model.
%
% USAGE:
%
%    FileString = generateFieldDescriptionFile(FileName)
%
% OPTIONAL INPUT:
%    FileName:      The FileName to write to. (default [CBTDIR filesep 'docs' filesep 'notes' filesep 'COBRAModelFields.md'])
%
% OUTPUT:
%    FileString:    The string written to the specified filename (or the default file)
%
% .. Author: - Thomas Pfau May 2017

global CBTDIR
if ~exist('FileName','var')
    if isempty(CBTDIR)
        initCobraToolbox(false); %Don't update the toolbox automatically
    end
    FileName = [CBTDIR filesep 'docs' filesep 'source' filesep 'notes' filesep 'COBRAModelFields.md'];
end

fieldProperties = getDefinedFieldProperties('Descriptions', true);

FileString = sprintf('# Fields in the model structure\n\n');
FileString = strjoin({FileString, sprintf(['Contents:\n',...
                                        '1. [Model Fields](#model-fields)\n',...
                    '2. [Field Support](#field-support)\n',...
					'3. [Model Specific Fields](#model-specific-fields)\n',...
					'4. [Annotation Definitions](#annotation-definitions)\n',...
					'\n',...
					'### Model Fields:\n',...
					'The following fields are defined in the COBRA toolbox. IF the field is present in a model, it should have the properties defined here and should be of the mentioned dimensions.\n '...
                    'The dimensions refer to m (the number of metabolites), n (the number of reactions), g (the number of genes) and c (the number of compartments).\n',...
					'\n',...
					'| Field Name | Dimension | Field Type | Field Description |\n',...
					'|---|---|---|---|\n'])},'');

for i=1:size(fieldProperties,1)
	string = sprintf('|`model.%s`| `%s` | %s | %s | \n',fieldProperties{i,1},fieldProperties{i,2},fieldProperties{i,3},fieldProperties{i,4});
	FileString = strjoin({FileString,string},'');
end
FileString = strjoin({FileString,sprintf(['### Model Specific Fields\n',...
					'Some models might contain additional model specific fields that are not defined COBRA model fields. These fields will commonly not be considered by COBRA toolbox methods, and using toolbox methods can render these fields inconsistent (e.g. if the number of reactions changes, a model specific field linked to reactions might have the wrong number of entries or the values might no longer correspond to the correct indices).\n',...
					'\n',...
					'### Field Support\n',...
					'All fields mentioned above are supported by COBRA Toolbox functions.'...
                    'Using COBRA Toolbox Functions will not make a model inconsistent, but manual modifications of fields might lead to an inconsistent model.\n',...
					'Use verifyModel(model) to determine, if the model is a valid COBRA Toolbox model.\n',...
					'\n',...
					'### Additional fields\n',...
					'Fields starting with met, rxn, comp, protein or gene that are not defined above, will be assumed to be annotation fields, and IO methods will try to map them to identifiers.org registered databases.'])},'');
fileID = fopen(FileName,'w');
fprintf(fileID,FileString);
fclose(fileID);


text:
function [ point, seed ] = p14_sample ( m, n, seed )

%*****************************************************************************80
%
%% P14_SAMPLE samples points from the region in problem 14.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real POINT(M,N), the coordinates of the points.
%
%    Output, integer SEED, a seed for the random number generator.
%
  box = [ ...
    100.0,  145.0; ...
   +634.0, +799.0 ]';

  reject = 0;

  for j = 1 : n

    while ( 1 )

      [ x, seed  ] = r8mat_uniform_01 ( m, 1, seed );
      x(1:m,1) = ( 1.0 - x(1:m,1) ) .* box(1:m,1) + x(1:m,1) .* box(1:m,2);

      inside = p14_inside ( m, 1, x );

      if ( inside )
        break
      end

      reject = reject + 1;

      if ( 30 * n + 10 <= reject )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'P14_SAMPLE - Fatal error!\n' );
        fprintf ( 1, '  Trying to generate point J = %d\n', j );
        fprintf ( 1, '  Number of rejections = %d\n', reject );
        fprintf ( 1, '  Rejection percentage = %f\n', ...
          ( 100 * reject ) / ( reject + j - 1 ) );
        r8vec_print ( m, x, '  Most recent rejected point: ' );
        error ( 'P14_SAMPLE - Fatal error!' );
      end

    end

    point(1:m,j) = x(1:m,1);

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Local Regression and Likelihood, Figure 2.7.
%
% Example gcvplot. First argument to gcvplot is a
% vector (or matrix) of smoothing parameters.
% Remaining arguments are passed directly to locfit().
%
% Author: Catherine Loader.
%
% NEED: cpplot.

load ethanol;

alp = (0.2:0.05:0.8)';
figure('Name','fig2_7: Example gcvplot');
gcvplot(alp,E,NOx);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%for subplots
function p=plots(x,y,z,img)
H1=abs(img);
colormap(hot)
subplot(x,y,z);
image(100*H1/max(max(H1)));
if z==2
    title('ISAR images using Harmonic Wavelets');
end
if z==4
    ylabel('Range')
end
if z==8
    xlabel('Doppler')
end


text:
function eventLog = ma_executeCoast_goto_tru(truTarget, initialState, eventNum, considerSoITransitions, soiSkipIds, refBody, massLoss, orbitDecay, celBodyData)
%ma_executeCoast_goto_tru Summary of this function goes here
%   Detailed explanation goes here
    bodyID = initialState(8);

    bodyInfo = getBodyInfoByNumber(bodyID, celBodyData);
    gmu = bodyInfo.gm;
    rVect = initialState(2:4)';
    vVect = initialState(5:7)';

    [sma, ecc, ~, ~, ~, truINI] = getKeplerFromState(rVect,vVect,gmu);
    meanMotion = computeMeanMotion(sma, gmu);

    if(isempty(refBody))
        refBody.id = -1;
    end

    if((abs(truTarget-truINI) < 1E-7 || abs(abs(truTarget-truINI)-2*pi) < 1E-7) && refBody.id==bodyInfo.id)
        eventLog = initialState;
        eventLog(:,13) = eventNum;
        return;
    end

    if(ecc < 1.0)
        if(abs(truTarget-2*pi) < 1E-7)
            truTarget = 2*pi;
        else
            truTarget = AngleZero2Pi(truTarget);
        end

        if(abs(truINI-2*pi) < 1E-7)
            if(truINI - 2*pi >= 0 || truINI==2*pi)
                truINI = 0;
            else
                truINI = 0;
            end
        else
            truINI = AngleZero2Pi(truINI);
        end
    else
        truINI = angleNegPiToPi(truINI);
        truTarget = angleNegPiToPi(truTarget);
    end

    if(abs(AngleZero2Pi(truTarget)-AngleZero2Pi(truINI)) <1E-6)
        truTarget = truINI;
    end

    if(ecc >= 1.0)
        parentBodyInfo = bodyInfo.getParBodyInfo(celBodyData);
        if(~isempty(parentBodyInfo))
            rSOI = getSOIRadius(bodyInfo, parentBodyInfo);
        else
            rSOI = Inf;
        end
        iniOrbit = [sma, ecc];
        iniHyTruMax = AngleZero2Pi(computeTrueAFromRadiusEcc(rSOI, iniOrbit(1), iniOrbit(2)));
        meanMotion = computeMeanMotion(iniOrbit(1), bodyInfo.gm);
        lbTA = -iniHyTruMax;
        ubTA = iniHyTruMax;

        bool1 = (ecc >= 1.0 && (truTarget)<(truINI));
        bool2 = (ecc >= 1.0 && truTarget>ubTA);
        bool3 = (ecc >= 1.0 && truTarget<lbTA);

        if(bool1 || bool2 || bool3)
            soITrans = findSoITransitions(initialState, Inf, soiSkipIds, massLoss, orbitDecay, celBodyData);
            if(~isempty(soITrans))
                SoITransEventLog = ma_executeCoast_goto_soi_trans(initialState, eventNum, Inf, soiSkipIds, massLoss, orbitDecay, celBodyData, soITrans);
                goToUTEventLog = ma_executeCoast_goto_tru(truTarget, SoITransEventLog(end,:), eventNum, true, soiSkipIds, refBody, massLoss, orbitDecay, celBodyData);
                eventLog = [SoITransEventLog; goToUTEventLog];
                return;
            else
                warnStr = ['Cannot go to target true anomaly (', num2str(rad2deg(truTarget)), ' deg).  Could not find an SoI transition, skipping.'];
                addToExecutionWarnings(warnStr, eventNum, bodyID, celBodyData);
                eventLog = initialState;
                eventLog(:,13) = eventNum;
                return;
            end
        end
    end

    if(considerSoITransitions)
        if(refBody.id == bodyID || refBody.id == -1)
            dt = getDtToTru(truTarget, sma, ecc, truINI, gmu);
            utTru = initialState(1) + dt;
        else
            dt = NaN;
            utTru = Inf;
        end

        if(not(isfinite(utTru)) || dt > 0)
            soITrans = findSoITransitions(initialState, utTru, soiSkipIds, massLoss, orbitDecay, celBodyData);
            SoITransEventLog = [];
            if(~isempty(soITrans) && min(soITrans(:,2)) < utTru)
                SoITransEventLog = ma_executeCoast_goto_soi_trans(initialState, eventNum, utTru, soiSkipIds, massLoss, orbitDecay, celBodyData, soITrans);
                goToUTEventLog = ma_executeCoast_goto_tru(truTarget, SoITransEventLog(end,:), eventNum, true, soiSkipIds, refBody, massLoss, orbitDecay, celBodyData);
            else
                goToUTEventLog = ma_executeCoast_goto_tru(truTarget, initialState, eventNum, false, soiSkipIds, refBody, massLoss, orbitDecay, celBodyData);
            end
            eventLog = [SoITransEventLog; goToUTEventLog];
        else
            eventLog = initialState;
        end

        return;
    end

    if(truTarget >= truINI)
        meanINI = computeMeanFromTrueAnom(truINI, ecc);
        if(abs(truTarget) < 1E-8)
            meanTarget = 0;
        elseif(abs(truTarget-2*pi) < 1E-8)
            meanTarget = 2*pi;
            if(meanINI < 0)
                meanINI = meanINI + 2*pi;
            end
        else
            meanTarget = computeMeanFromTrueAnom(truTarget, ecc);
        end

        dM = meanTarget-meanINI;
        dt = dM/meanMotion;

        eventLog = ma_executeCoast_goto_dt(dt, initialState, eventNum, considerSoITransitions, soiSkipIds, massLoss, orbitDecay, celBodyData);
    else
        eventLog1 = ma_executeCoast_goto_tru(2*pi, initialState, eventNum, considerSoITransitions, soiSkipIds, refBody, massLoss, orbitDecay, celBodyData);

        bodyID = eventLog1(end,8);
        bodyInfo = getBodyInfoByNumber(bodyID, celBodyData);
        gmu = bodyInfo.gm;
        [sma, ecc, inc, raan, arg, ~] = getKeplerFromState(eventLog1(end,2:4),eventLog1(end,5:7),gmu);

        if(truTarget>=1E-6)
            newTruINI = 1E-6;
            [rVect,vVect]=getStatefromKepler(sma, ecc, inc, raan, arg, newTruINI, gmu); %set true anomaly to 0 to avoid having issues w/ infinite recursion
            eventLog1(end,2:4) = reshape(rVect,1,3);
            eventLog1(end,5:7) = reshape(vVect,1,3);

            eventLog2 = ma_executeCoast_goto_tru(truTarget, eventLog1(end,:), eventNum, considerSoITransitions, soiSkipIds, refBody, massLoss, orbitDecay, celBodyData);
        else
            eventLog2 = [];
        end

        eventLog = [eventLog1; eventLog2];
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


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

% DTFT properties

% Differentiation in frequency

syms n
x=0.8.^n
L=symsum(n/j*x*exp(-j*w*n),n,0,inf);
subplot(211);
w1=-pi:.01:pi;
Li=subs(L,w,w1);
plot(w1,abs(Li));
legend('magnitude of DTFT of left part')
subplot(212);
plot(w1,angle(Li))
legend('angle of DTFT of left part')

figure
syms n w
X=1/(1-0.8*exp(-j*w));
R=diff(X,w);
subplot(211);
w1=-pi:.01:pi;
Ri=subs(R,w,w1);
plot(w1,abs(Ri));
legend('magnitude of right part')
subplot(212);
plot(w1,angle(Ri))
legend('angle of right part')


text:
function b = tril(a,k)
% function B = tril(A,K)
%
% DESCRIPTION
%   Extract lower triangular part of a polynomial matrix
%
% INPUTS
%   A: polynomial
%   K: integer
%
% OUTPUTS
%   B: polynomial
%
% SYNTAX
%   B=tril(A,K)
%     B contains the elements on and below the K-th diagonal of A. K=0
%     is the main diagonal, K>0 is above the main diagonal, and K<0 is
%     below the main diagonal.
%   B=tril(A)
%     This is the same as B=tril(A,0).
%
% See also triu, diag

% 11/24/2009: PJS  Initial Coding

if nargin==1
    k=0;
end

sza = size(a);
% a is a matrix --> b is a vector
if k>=sza(2)
    b = a;
    return
elseif k <= -sza(1)
    b = polynomial(zeros(sza));
    return;
end

% Use double/tril to find indices to zero out
idx = find( ~tril(ones(sza),k) );

% Assign values
b = a;
L.type = '()';
L.subs = {idx};
b = subsasgn(b,L,0);

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
classdef MrCopyData < handle
    % Provides a common clone-method for all object classes, plus
    % generalized find, compare and print-capabilities
    %
    % based on a posting by Volkmar Glauche
    % http://www.mathworks.com/matlabcentral/fileexchange/22965-clone-handle-object-using-matlab-oop
    %
    % heavily modified and extended by Lars Kasper and Saskia Klein

    % Author:   Saskia Klein & Lars Kasper
    % Created:  2010-04-15
    % Copyright (C) 2014 Institute for Biomedical Engineering
    %                    University of Zurich and ETH Zurich
    %
    % This file is part of the TAPAS UniQC Toolbox, which is released
    % under the terms of the GNU General Public Licence (GPL), version 3.
    % You can redistribute it and/or modify it under the terms of the GPL
    % (either version 3 or, at your option, any later version).
    % For further details, see the file COPYING or
    %  <http://www.gnu.org/licenses/>.
    %

    properties
    end
    methods
        function [this, unusedArgs] = MrCopyData(varargin)
            % Constructor for MrCopyData
            % either
            %
            %   this = MrCopyData for object with default values
            %
            %   OR
            %
            %   this = MrCopyData('empty') to create an object with all
            %   values set to [];
            %
            %   OR
            %   this = MrCopyData('param_name1', param_value1, 'param_name2', param_value2 )
            %          set of parameter names and values given, e.g.
            %          MrCopyData('dyn', 1)
            unusedArgs = {};
            if nargin
                if strcmpi(varargin{1}, 'empty')
                    this.clear();
                else
                    for cnt = 1:nargin/2 % save them to object properties
                        if isprop(this, varargin{2*cnt-1})
                        this.(varargin{2*cnt-1}) = varargin{2*cnt};
                        else % return unused props in a struct
                            unusedArgs = [unusedArgs, varargin(2*cnt-1:2*cnt)];
                        end
                    end
                end
            end
        end
    end
end

author: translationalneuromodeling

repo: tapas


text:
function [inface, outface]=innersurf(node,face,outface)
%
% outface=innersurf(node,face,outface)
%
% extract the interior triangles (shared by two enclosed compartments) of a complex surface
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%    node:  node coordinates
%    face:  surface triangle list
%    outface: (optional) the exterior triangle list, if not given, will
%           be computed using outersurf().
%
% output:
%    inface: the collection of interior triangles of the surface mesh
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if(nargin<3)
  outface=outersurf(node,face);
end

[I,J]=ismember(sort(face,2),sort(outface,2),'rows');

inface=face(I==0,:);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function test12
%TEST12 test cs_qr and compare with svd
%
% Example:
%   test12
% See also: testall

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

fprintf ('test 12\n') ;
rand ('state',0) ;
% A = rand (3,4)

for trial = 1:100
    m = fix (100 * rand (1)) ;
    n = fix (100 * rand (1)) ;
    d = .1 * rand (1) ;
    A = sprandn (m,n,d) ;
    fprintf ('m %d n %d nnz %d\n', m, n, nnz(A)) ;
    if (m < n)
        continue ;
    end
    if (m == 0 | n == 0)                                                    %#ok
        continue ;
    end
    % save A A
    fprintf ('[ ') ;
    [V,Beta,p,R] = cs_qr (A) ;
    % [Q,R] = svd (full(A)) ;
    fprintf (']\n') ;

    s1 = svd (full (A)) ;
    s2 = svd (full (R)) ;
    s2 = s2 (1:length(s1)) ;
    err = norm (s1-s2) ;
    if (length (s1) > 1)
        err = err / s1 (1) ;
    end
    fprintf ('err %g\n', err) ;
    if (err > 1e-12)
        error ('!') ;
    end
end