text:
function [esTSNR,esHRNR]=WienerNoiseReduction(ns,fs,IS)

% [esTSNR,esHRNR]=WIENERNOISEREDUCTION(ns,fs,IS)
%
% Title  :      Wiener Noise Suppressor with TSNR & HRNR algorithms
%
% Description : Wiener filter based on tracking a priori SNR using Decision-Directed
%               method, proposed by Plapous et al 2006. The two-step noise reduction
%               (TSNR) technique removes the annoying reverberation effect while
%               maintaining the benefits of the decision-directed approach. However,
%               classic short-time noise reduction techniques, including TSNR, introduce
%               harmonic distortion in the enhanced speech. To overcome this problem, a
%               method called harmonic regeneration noise reduction (HRNR)is implemented
%               in order to refine the a priori SNR used to compute a spectral gain able
%               to preserve the speech harmonics.
%
%
% Reference :   Plapous, C.; Marro, C.; Scalart, P., "Improved Signal-to-Noise Ratio
%               Estimation for Speech Enhancement", IEEE Transactions on Audio, Speech,
%               and Language Processing, Vol. 14, Issue 6, pp. 2098 - 2108, Nov. 2006
%
% Input Parameters :
%   ns          Noisy speech
%   fs          Sampling frequency (in Hz)
%   IS          Initial Silence (or non-speech activity) Period (in number of samples)
%
% Output Parameters : enhanced speech
%   esTSNR      enhanced speech with the Two-Step Noise Reduction method
%   esHNRN      enhanced speech with the Harmonic Regeneration Noise Reduction method
%
%Author :       LIU Ming, 2008
%Modified :     SCALART Pascal october, 2008
%
%

%% ------- input noisy speech  --------

l = length(ns);
s=ns;

wl = fix(0.020*fs)    % window length is 20 ms
NFFT=2*wl             % FFT size is twice the window length
hanwin = hanning(wl);

if (nargin<3 | isstruct(IS))
    IS=10*wl;             %Initial Silence or Noise Only part in samples (= ten frames)
end
%% -------- compute noise statistics ----------

nsum = zeros(NFFT,1);

count = 0;
    for m = 0:IS-wl
     nwin = s(m+1:m+wl).*hanwin;
      nsum = nsum + abs(fft(nwin,NFFT)).^2;
     count = count + 1;
    end

d= (nsum)/count;

%% --------- main algorithm ---------------
SP = 0.25;      % Shift percentage is 50 % Overlap-Add method works good with this value
normFactor=1/SP;
overlap = fix((1-SP)*wl); % overlap between sucessive frames
offset = wl - overlap;
max_m = fix((l-NFFT)/offset);

zvector = zeros(NFFT,1);
oldmag = zeros(NFFT,1);
news = zeros(l,1);

phasea=zeros(NFFT,max_m);
xmaga=zeros(NFFT,max_m);
tsnra=zeros(NFFT,max_m);
newmags=zeros(NFFT,max_m);

alpha = 0.98;

%Iteration to remove noise

%% --------------- TSNR ---------------------
for m = 0:max_m
   begin = m*offset+1;
   iend = m*offset+wl;
   speech = ns(begin:iend);       %extract speech segment
   winy = hanwin.*speech;   %perform hanning window
   ffty = fft(winy,NFFT);          %perform fast fourier transform
   phasey = angle(ffty);         %extract phase
   phasea(:,m+1)=phasey;       %for HRNR use
   magy = abs(ffty);             %extract magnitude
   xmaga(:,m+1)= magy;           %for HRNR use
   postsnr = ((magy.^2) ./ d)-1 ;      %calculate a posteriori SNR
   postsnr=max(postsnr,0.1);  % limitation to prevent distorsion

   %calculate a priori SNR using decision directed approach
   eta = alpha * ( (oldmag.^2)./d ) + (1-alpha) * postsnr;
   newmag = (eta./(eta+1)).*  magy;

   %calculate TSNR
   tsnr = (newmag.^2) ./ d;
   Gtsnr = tsnr ./ (tsnr+1);         %gain of TSNR
   tsnra(:,m+1)=Gtsnr;
   %Gtsnr=max(Gtsnr,0.1);
   Gtsnr = gaincontrol(Gtsnr,NFFT/2);

      %for HRNR use
   newmag = Gtsnr .* magy;
   newmags(:,m+1) = newmag;     %for HRNR use
   ffty = newmag.*exp(i*phasey);
   oldmag = abs(newmag);
   news(begin:begin+NFFT-1) = news(begin:begin+NFFT-1) + real(ifft(ffty,NFFT))/normFactor;
end

esTSNR=news;
%% --------------- HRNR -----------------------

%non linearity
newharm= max(esTSNR,0);
news = zeros(l,1);
%
for m = 0:max_m
   begin = m*offset+1;
   iend = m*offset+wl;

   nharm = hanwin.*newharm(begin:iend);
   ffth = abs(fft(nharm,NFFT));          %perform fast fourier transform

   snrham= ( (tsnra(:,m+1)).*(abs(newmags(:,m+1)).^2) + (1-(tsnra(:,m+1))) .* (ffth.^2) ) ./d;

   newgain= (snrham./(snrham+1));
   %newgain=max(newgain,0.1);

   newgain=gaincontrol(newgain,NFFT/2);

   newmag = newgain .*  xmaga(:,m+1);

   ffty = newmag.*exp(i*phasea(:,m+1));

   news(begin:begin+NFFT-1) = news(begin:begin+NFFT-1) + real(ifft(ffty,NFFT))/normFactor;
end;
%Output
esHRNR=news;

figure;
[B,f,T] = specgram(ns,NFFT,fs,hanning(wl),wl-10);
imagesc(T,f,20*log10(abs(B)));axis xy;colorbar
title(['Spectrogram - noisy speech'])
xlabel('Time (sec)');ylabel('Frequency (Hz)');

figure;
[B,f,T] = specgram(esTSNR,NFFT,fs,hanning(wl),wl-10);
imagesc(T,f,20*log10(abs(B)));axis xy;colorbar
title(['Spectrogram - output speech TSNR'])
xlabel('Time (sec)');ylabel('Frequency (Hz)');

figure;
[B,f,T] = specgram(esHRNR,NFFT,fs,hanning(wl),wl-10);
imagesc(T,f,20*log10(abs(B)));axis xy;colorbar
title(['Spectrogram - output speech HRNR'])
xlabel('Time (sec)');ylabel('Frequency (Hz)');

function        NewGain=gaincontrol(Gain,ConstraintInLength)
%
%Title  : Additional Constraint on the impulse response
%         to ensure linear convolution property
%
%
%Description :
%
% 1- The time-duration of noisy speech frame is equal to L1 samples.
%
% 2- This frame is then converted in the frequency domain
%       by applying a short-time Fourier transform of size NFFT leading
%       to X(wk) k=0,...,NFFT-1 when NFFT is the FFT size.
%
% 3- The estimated noise reduction filter is G(wk) k=0,1,...,NFFT-1
%       leading to an equivalent impulse response g(n)=IFFT[G(wk)]
%       of length L2=NFFT
%
% 4- When applying the noise reduction filter G(wk) to the noisy
%       speech spectrum X(wk), the multiplication S(wk)=G(wk)X(wk) is
%       equivalent to a convolution in the time domain. So the
%       time-duration of the enhanced speech s(n) should be equal to
%       Ltot=L1+L2-1.
%
% 5- If the length Ltot is greater than the time-duration of the IFFT[S(wk)]
%       the a time-aliasing effect will appear.
%
% 6- To overcome this phenomenon, the time-duration L2 of the equivalent
%       impulse response g(n) should be chosen such that Ltot = L1 + L2 -1 <= NFFT
%       => L2 <= NFFT+1-Ll
%
%       here we have NFFT=2*Ll so we should have L2 <= Ll+1. I have made
%       the following choice : the time-duration of g(n) is limited to
%       L2=NFFT/2=L1 (see lines 88 and 192)
%
%Author : SCALART Pascal
%
%October  2008
%

meanGain=mean(Gain.^2);
NFFT=length(Gain);

L2=ConstraintInLength;

win=hamming(L2);

% Frequency -> Time
% computation of the non-constrained impulse response
ImpulseR=real(ifft(Gain));

% application of the constraint in the time domain
ImpulseR2=[ImpulseR(1:L2/2).*win(1+L2/2:L2);zeros(NFFT-L2,1);ImpulseR(NFFT-L2/2+1:NFFT).*win(1:L2/2)];

% Time -> Frequency
NewGain=abs(fft(ImpulseR2,NFFT));

meanNewGain=mean(NewGain.^2);

NewGain=NewGain*sqrt(meanGain/meanNewGain); % normalisation to keep the same energy (if white r.v.)

author: jtkim-kaist

repo: Speech-enhancement


text:
% function inlap
% inlap.m
%
% make dialog interface for seislap
%
% Last change 7/95               A.Allmann

%
report_this_filefun(mfilename('fullpath'));

global ptime
global seis_inpu1 seis_inpu2 seis_inpu3 qui1
global repeat_button lapf maepi

%set catalog
if size(newcat)==0
    newcat=a;
end

%initial values
figure_w_normalized_uicontrolunits(mess)
clf
set(gca,'visible','off')
set(gcf,'Units','pixel','NumberTitle','off','Name','Input Parameters for Seislap');

set(gcf,'pos',[ wex  wey welx+100 wely+50])

bev=find(newcat.Magnitude==max(newcat.Magnitude)); %biggest events in catalog

%default values of input parameters
ldx=100;
latt=newcat(bev(1),2);
longt=newcat(bev(1),1);
Mmin=3;
ldepth=newcat(bev(1),7);
if quie==1
    qui1=1;
    tlap=100;
    binlength=10;
elseif quie==2
    qui1=2;
    tlap=10;
    binlength=1;
end

% creates a dialog box to input some parameters
%
txt1 = text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0. 0.96 0 ],...
    'Rotation',0 ,...
    'FontWeight','bold',...
    'FontSize',ZmapGlobal.Data.fontsz.m ,...
    'String','Minimum Magnitude: ');

txt2 = text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0. 0.81 0 ],...
    'Rotation',0 ,...
    'FontWeight','bold',...
    'FontSize',ZmapGlobal.Data.fontsz.m ,...
    'String','Please input DX in km ');

txt3 = text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0. 0.66 0 ],...
    'Rotation',0 ,...
    'FontWeight','bold',...
    'FontSize',ZmapGlobal.Data.fontsz.m ,...
    'String','Please input overlap time in days:');

txt4 = text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0. 0.51 0 ],...
    'Rotation',0 ,...
    'FontWeight','bold',...
    'FontSize',ZmapGlobal.Data.fontsz.m ,...
    'String','Longitude of point of interest:');
txt5 = text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0. 0.36 0 ],...
    'Rotation',0 ,...
    'FontWeight','bold',...
    'FontSize',ZmapGlobal.Data.fontsz.m ,...
    'String','Lattitude of point of interest:');

txt6 = text(...
    'Color',[0 0 0 ],...
    'EraseMode','normal',...
    'Position',[0. 0.21 0 ],...
    'Rotation',0 ,...
    'FontWeight','bold',...
    'FontSize',ZmapGlobal.Data.fontsz.m ,...
    'String','Binlength in days:');

inp1_field  = uicontrol('Style','edit',...
    'Position',[.85 .87 .13 .08],...
    'Units','normalized','String',num2str(Mmin),...
    'Callback','Mmin=str2double(get(inp1_field,''String''));set(inp1_field,''String'',num2str(Mmin));');

inp2_field  = uicontrol('Style','edit',...
    'Position',[.85 .75 .13 .08],...
    'Units','normalized','String',num2str(ldx),...
    'Callback','ldx=str2double(get(inp2_field,''String''));set(inp2_field,''String'',num2str(ldx));');

inp3_field=uicontrol('Style','edit',...
    'Position',[.85 .63 .13 .08],...
    'Units','normalized','String',num2str(tlap),...
    'Callback','tlap=str2double(get(inp3_field,''String''));     set(inp3_field,''String'',num2str(tlap));');

inp4_field=uicontrol('Style','edit',...
    'Position',[.85 .51 .13 .08],...
    'Units','normalized','String',num2str(longt),...
    'Callback','longt=str2double(get(inp4_field,''String'')); set(inp4_field,''String'',num2str(longt));');

inp5_field=uicontrol('Style','edit',...
    'Position',[.85 .39 .13 .08],...
    'Units','normalized','String',num2str(latt),...
    'Callback','latt=str2double(get(inp5_field,''String'')); set(inp5_field,''String'',num2str(latt));');

inp6_field=uicontrol('Style','edit',...
    'Position',[.85 .27 .13 .08],...
    'Units','normalized','String',num2str(binlength),...
    'Callback','binlength=str2double(get(inp6_field,''String'')); set(inp6_field,''String'',num2str(binlength));');

close_button=uicontrol('Style','Pushbutton',...
    'Position', [.7 .05 .15 .12 ],...
    'Units','normalized','Callback','welcome','String','Cancel');

inflap_button=uicontrol('Style','Pushbutton',...
    'Position', [.47 .05 .15 .12 ],...
    'Units','normalized','Callback','clinfo(13)','String','Info');

go_button=uicontrol('Style','Pushbutton',...
    'Position',[.25 .05 .15 .12 ],...
    'Units','normalized',...
    'Callback','tlap=str2num(get(inp3_field,''String''));ldx=str2num(get(inp2_field,''String''));binlength=str2num(get(inp6_field,''String''));latt=str2num(get(inp5_field,''String''));longt=str2num(get(inp4_field,''String''));Mmin=str2num(get(inp1_field,''String''));clf;welcome;if quie==1 lap1=seislap2(ldx,tlap,binlength,longt,latt,Mmin,ldepth,1);elseif quie==2 lap1=seislap2(ldx,tlap,binlength,longt,latt,Mmin,ldepth,2);end;',...
    'String','Go');

%freeze_button = uicontrol(...
%'BackgroundColor',[ 0.7 0.7 0.7 ],...
%'Callback','fre = get(freeze_button,''Value'');',...
%'ForegroundColor',[ 0 0 0 ],...
%'Position',[ 0.25 0.30 0.4 0.15 ],...
%'String','Freeze Colorbar? ',...
%'Style','checkbox',...
%'Units','normalized',...
%'Visible','on');
%

set(gcf,'visible','on');watchoff

author: CelsoReyes


text:
function det = r8pbu_det ( n, mu, a_lu )

%*****************************************************************************80
%
%% R8PBU_DET computes the determinant of a matrix factored by R8PBU_FA.
%
%  Discussion:
%
%    The R8PBU storage format is for a symmetric positive definite band matrix.
%
%    To save storage, only the diagonal and upper triangle of A is stored,
%    in a compact diagonal format that preserves columns.
%
%    The diagonal is stored in row MU+1 of the array.
%    The first superdiagonal in row MU, columns 2 through N.
%    The second superdiagonal in row MU-1, columns 3 through N.
%    The MU-th superdiagonal in row 1, columns MU+1 through N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 March 2004
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Bunch, Moler, Stewart,
%    LINPACK User's Guide,
%    SIAM, Philadelphia, 1979.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, integer MU, the number of superdiagonals of the matrix.
%    MU must be at least 0 and no more than N-1.
%
%    Input, real A_LU(MU+1,N), the LU factors from R8PBU_FA.
%
%    Output, real DET, the determinant of the matrix.
%
  det = prod ( a_lu(mu+1,1:n) );
  det = det * det;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
clc;clear;
close all;
Angle = [90,180];                               %Angle,Options:90,135,180,270 plus:0,45,225,315
if(exist('SetupStruc','var'))
    [s,sOri,Unmix_s,~] = readData(Angle,ISM_setup);
else
    [s,sOri,Unmix_s,SetupStruc] = readData(Angle,ISM_setup);    % 's' is the muli-channel mixture signal,size:Length*7,
end                                            % 'sOri' is the each single singal of the 1st channel corresponding to the angle,size:Length*N
method = {'DSB'                1;
          'DSB_Mask'           0;
          'MVDR'               0;
          'MVDR_ESB'           0;
          'MVDR_AESB'          0;
          'MVDR_PCA'           0;
          'MVDR_Search'        0;
          'MVDR_AESB_Search'   0;
          'LCMV'               0;
          'LCMV_ESB'           0;
          'LCMV_AESB'          0;
          'LCMV_Search'        0;
          'ICA_funda'          0;
          'ICA_initial'        0;
          'ICA_Sawada'         0;
          'FastICA_HO_Sawada'  0;
          'IVA'                0;
          'IVA_woDR'           0;
          'AuxIVA'             0;
          'OverIVA'            0;
          'OverILRMA'          0;
          'ILRMA_woDR'         0;
          'ILRMA'              0;
          'ILRMA_PF'           0;
          'FastMNMF1'          0;
          'FastMNMF2'          0;
          'FastFCA_AS'         0; %%%% remaining to be written
          'maxSNR'             0;
          %%%%%%% Dereverberation
          'WPE'                0;
          %%%%%%% Compound method
          'cGMM_maxSNR'        0;
          'ALL'                0
          };
SetupStruc.unS = Unmix_s;
[Re,SetupStruc] = Process(method,s,SetupStruc);
[Me] = Cal_metrics(Re,Unmix_s,sOri,SetupStruc);

repo: Voice-Separation-and-Enhancement


text:
function u = nvecs(X,n,r,opts)
%NVECS Compute the leading mode-n vectors for a tensor.
%
%   U = NVECS(X,n,r) computes the r leading eigenvalues of Xn*Xn'
%   (where Xn is the mode-n matricization of X), which provides
%   information about the mode-n fibers. In two-dimensions, the r
%   leading mode-1 vectors are the same as the r left singular vectors
%   and the r leading mode-2 vectors are the same as the r right
%   singular vectors.
%
%   U = NVECS(X,n,r,OPTS) specifies options:
%   OPTS.eigsopts: options passed to the EIGS routine [struct('disp',0)]
%   OPTS.flipsign: make each column's largest element positive [true]
%
%   See also TENSOR, TENMAT, EIGS.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if ~exist('opts','var')
    opts = struct;
end

if isfield(opts,'eigsopts')
    eigsopts = opts.eigsopts;
else
    eigsopts.disp = 0;
end

Xn = double(tenmat(X,n));
Y = Xn*Xn';
warning off;
[u,d] = eigs(Y, r, 'LM', eigsopts);
warning on;
if isfield(opts,'flipsign')
    flipsign = opts.flipsign;
else
    flipsign = true;
end

if flipsign
    % Make the largest magnitude element be positive
    [val,loc] = max(abs(u));
    for i = 1:r
        if u(loc(i),i) < 0
            u(:,i) = u(:,i) * -1;
        end
    end
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function cost = costSymbolSeparation(par, sampleRate, iq)

symbols = computeSymbolsFromParams(par./[1 1 100 100], sampleRate, iq);

%%
aa = mod(angle(symbols),(pi/4));
bb = mod(abs(symbols),.25);

%%
aa = min([aa; pi/4-aa]);
bb = min([bb; .25-bb]);

%%
cost = sum([aa.^2, bb.^2]);

author: Sable

repo: mcbench-benchmarks


text:
function varargout = lvd_EditActionSetTankToTankConnFlowRateGUI(varargin)
% LVD_EDITACTIONSETTANKTOTANKCONNFLOWRATEGUI MATLAB code for lvd_EditActionSetTankToTankConnFlowRateGUI.fig
%      LVD_EDITACTIONSETTANKTOTANKCONNFLOWRATEGUI, by itself, creates a new LVD_EDITACTIONSETTANKTOTANKCONNFLOWRATEGUI or raises the existing
%      singleton*.
%
%      H = LVD_EDITACTIONSETTANKTOTANKCONNFLOWRATEGUI returns the handle to a new LVD_EDITACTIONSETTANKTOTANKCONNFLOWRATEGUI or the handle to
%      the existing singleton*.
%
%      LVD_EDITACTIONSETTANKTOTANKCONNFLOWRATEGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in LVD_EDITACTIONSETTANKTOTANKCONNFLOWRATEGUI.M with the given input arguments.
%
%      LVD_EDITACTIONSETTANKTOTANKCONNFLOWRATEGUI('Property','Value',...) creates a new LVD_EDITACTIONSETTANKTOTANKCONNFLOWRATEGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before lvd_EditActionSetTankToTankConnFlowRateGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to lvd_EditActionSetTankToTankConnFlowRateGUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help lvd_EditActionSetTankToTankConnFlowRateGUI

% Last Modified by GUIDE v2.5 21-Jan-2019 13:19:19

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @lvd_EditActionSetTankToTankConnFlowRateGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @lvd_EditActionSetTankToTankConnFlowRateGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before lvd_EditActionSetTankToTankConnFlowRateGUI is made visible.
function lvd_EditActionSetTankToTankConnFlowRateGUI_OpeningFcn(hObject, eventdata, handles, varargin)
    % This function has no output args, see OutputFcn.
    % hObject    handle to figure
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)
    % varargin   command line arguments to lvd_EditActionSetTankToTankConnFlowRateGUI (see VARARGIN)

    % Choose default command line output for lvd_EditActionSetTankToTankConnFlowRateGUI
    handles.output = hObject;

    centerUIFigure(hObject);

    action = varargin{1};
    setappdata(hObject,'action',action);

    lv = varargin{2};
    setappdata(hObject,'lv',lv);

    populateGUI(handles, action, lv);

    % Update handles structure
    guidata(hObject, handles);

    % UIWAIT makes lvd_EditActionSetTankToTankConnFlowRateGUI wait for user response (see UIRESUME)
    uiwait(handles.lvd_EditActionSetTankToTankConnFlowRateGUI);

function populateGUI(handles, action, lv)
    [t2TConnStr, t2TConns] = lv.getTankToTankConnectionsListBoxStr();
    set(handles.refConnCombo,'String',t2TConnStr);

    if(not(isempty(action.conn)))
        ind = find(t2TConns == action.conn,1,'first');
    else
        ind = [];
    end

    if(not(isempty(ind)))
        set(handles.refConnCombo,'Value',ind);
    end

    handles.flowRateText.String = fullAccNum2Str(action.flowRateToSet);

% --- Outputs from this function are returned to the command line.
function varargout = lvd_EditActionSetTankToTankConnFlowRateGUI_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
    if(isempty(handles))
        varargout{1} = false;
    else
        action = getappdata(hObject,'action');
        lv = getappdata(hObject,'lv');

        [~, t2TConns] = lv.getTankToTankConnectionsListBoxStr();
        ind = get(handles.refConnCombo,'Value');
        action.conn = t2TConns(ind);

        action.flowRateToSet = str2double(handles.flowRateText.String);

        varargout{1} = true;
        close(handles.lvd_EditActionSetTankToTankConnFlowRateGUI);
    end

% --- Executes on button press in saveAndCloseButton.
function saveAndCloseButton_Callback(hObject, eventdata, handles)
% hObject    handle to saveAndCloseButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    errMsg = validateInputs(handles);

    if(isempty(errMsg))
        uiresume(handles.lvd_EditActionSetTankToTankConnFlowRateGUI);
    else
        msgbox(errMsg,'Invalid Inputs','error');
    end

function errMsg = validateInputs(handles)
    errMsg = {};

    val = str2double(get(handles.flowRateText,'String'));
    enteredStr = get(handles.flowRateText,'String');
    numberName = 'Flow Rate';
    lb = 0;
    ub = Inf;
    isInt = false;
    errMsg = validateNumber(val, numberName, lb, ub, isInt, errMsg, enteredStr);

% --- Executes on button press in cancelButton.
function cancelButton_Callback(hObject, eventdata, handles)
% hObject    handle to cancelButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    close(handles.lvd_EditActionSetTankToTankConnFlowRateGUI);

% --- Executes on selection change in refConnCombo.
function refConnCombo_Callback(hObject, eventdata, handles)
% hObject    handle to refConnCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns refConnCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from refConnCombo

% --- Executes during object creation, after setting all properties.
function refConnCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to refConnCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in stateCombo.
function stateCombo_Callback(hObject, eventdata, handles)
% hObject    handle to stateCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns stateCombo contents as cell array
%        contents{get(hObject,'Value')} returns selected item from stateCombo

% --- Executes during object creation, after setting all properties.
function stateCombo_CreateFcn(hObject, eventdata, handles)
% hObject    handle to stateCombo (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: popupmenu controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on key press with focus on lvd_EditActionSetTankToTankConnFlowRateGUI or any of its controls.
function lvd_EditActionSetTankToTankConnFlowRateGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to lvd_EditActionSetTankToTankConnFlowRateGUI (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'return'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'enter'
            saveAndCloseButton_Callback(handles.saveAndCloseButton, [], handles);
        case 'escape'
            close(handles.lvd_EditActionSetTankToTankConnFlowRateGUI);
    end

function flowRateText_Callback(hObject, eventdata, handles)
% hObject    handle to flowRateText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of flowRateText as text
%        str2double(get(hObject,'String')) returns contents of flowRateText as a double
    newInput = get(hObject,'String');
    newInput = attemptStrEval(newInput);
    set(hObject,'String', newInput);

% --- Executes during object creation, after setting all properties.
function flowRateText_CreateFcn(hObject, eventdata, handles)
% hObject    handle to flowRateText (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [net, opts] = get_model(opts, addFC)
if nargin < 2, addFC = true; end

t0 = tic;
modelFunc = str2func(sprintf('models.%s', opts.arch));
[net, opts, in_name, in_dim] = modelFunc(opts);
logInfo('%s in %.2fs', opts.arch, toc(t0));

% + FC layer
if addFC
    convobj = dagnn.Conv('size', [1 1 in_dim opts.dim], ...
        'pad', 0, 'stride', 1, 'hasBias', true);
    params = convobj.initParams();
    net.addLayer('fc', convobj, {in_name}, {'logits'}, {'fc_w', 'fc_b'});
    p1 = net.getParamIndex('fc_w');
    p2 = net.getParamIndex('fc_b');
    net.params(p1).value = params{1};
    net.params(p2).value = params{2};
    net.params(p1).learningRate = opts.lrmult;
    net.params(p2).learningRate = opts.lrmult;

    in_name = 'logits';
    in_dim  = opts.dim;
end

% + l2 normalization layer
net.addLayer('L2norm', dagnn.LRN('param', [2*in_dim, 0, 1, 0.5]), ...
    {in_name}, {'feats_l2'});

% + loss layer
lossobj = str2func(opts.obj);
net.addLayer('loss', lossobj('opt', opts), {'feats_l2', 'labels'}, {'objective'});

% print
if 0
    net.print({'data', [opts.imageSize opts.imageSize 3 opts.batchSize]}, ...
        'MaxNumColumns', 4, 'Layers', [], 'Parameters', []);
end

end

author: kunhe

repo: FastAP-metric-learning


text:
% MORLET_2D_NODC computes the 2-D elliptic Morlet filter given a set of
%    parameters in Fourier domain
%
% Usage
%    gab = MORLET_2D_NODC(N, M, sigma, slant, xi, theta, offset)
%
% Input
%    N (numeric): Width of the filter.
%    M (numeric): Height of the filter.
%    sigma (numeric): Standard deviation of the envelope.
%    slant (numeric): Eccentricity of the elliptic envelope.
%       (the smaller slant, the larger angular resolution).
%    xi (numeric): The frequency peak.
%    theta (numeric): Orientation in radians of the filter.
%    offset (numeric, optional): 2-D vector reprensting the offset location
%       (default [0 0]).
%
% Output
%    gab (numeric): N-by-M matrix representing the gabor filter in spatial
%       domain.
%
% Description
%    Compute a Morlet wavelet in Fourier domain.
%
%    Morlet wavelets have a 0 DC component.
%
% See also
%    GABOR_2D, MORLET_2D_PYRAMID

function gab = morlet_2d_noDC(N, M, sigma, slant, xi, theta, offset)

	if ~exist('offset','var')
		offset = [0, 0];
	end
	[x , y] = meshgrid(1:M, 1:N);

	x = x - ceil(M/2) - 1 - offset(1);
	y = y - ceil(N/2) - 1 - offset(2);

	Rth = rotation_matrix_2d(theta);
	A = Rth\ [1/sigma^2, 0 ; 0 slant^2/sigma^2] * Rth ;
	s = x.* ( A(1,1)*x + A(1,2)*y) + y.*(A(2,1)*x + A(2,2)*y ) ;

	%normalize sucht that the maximum of fourier modulus is 1

	gaussian_envelope = exp( - s/2);
	oscilating_part = gaussian_envelope .* exp(1i*(x*xi*cos(theta) + y*xi*sin(theta)));
	K = sum(oscilating_part(:)) ./ sum(gaussian_envelope(:));
	gabc = oscilating_part - K.*gaussian_envelope;

	gab=1/(2*pi*sigma^2/slant^2)*fftshift(gabc);

end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function ft_write_event(filename, event, varargin)

% FT_WRITE_EVENT writes an event structure to a file, a message daemon
% listening on a network socked, or to another computer connected through
% the serial port.
%
% Use as
%   ft_write_event(filename, event, ...)
%
% The first argument is a string containing the filename. The second
% argument is a structure with the event. Multiple events can be
% represented as a structure array.
%
% Events are represented as
%   event.type      string
%   event.sample    expressed in samples, the first sample of a recording is 1
%   event.value     number or string
%   event.offset    expressed in samples
%   event.duration  expressed in samples
%   event.timestamp expressed in timestamp units, which vary over systems (optional)
%
% Events can also be written to special communication streams
% by specifying the target as URI instead of a filename. Supported are
%   buffer://<host>:<port>
%   fifo://<filename>
%   tcp://<host>:<port>
%   udp://<host>:<port>
%   mysql://<user>:<password>@<host>:<port>
%   rfb://<password>@<host>:<port>
%   serial:<port>?key1=value1&key2=value2&...
%   rfb://<password>@<host>:<port>
%
% See also FT_READ_HEADER, FT_READ_DATA, FT_READ_EVENT, FT_WRITE_DATA

% Copyright (C) 2007-2010 Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

global event_queue   % for fcdc_global
global db_blob       % for fcdc_mysql
if isempty(db_blob)
  db_blob = 0;
end

if iscell(filename)
  % use recursion to write to multiple event targets
  for i=1:numel(filename)
    ft_write_event(filename{i}, event, varargin);
  end
  return
end

% set the defaults
eventformat = ft_getopt(varargin, 'eventformat', ft_filetype(filename));
swapping    = ft_getopt(varargin, 'swapping', 'native');
append      = ft_getopt(varargin, 'append', 'yes');
maxqlength  = ft_getopt(varargin, 'maxqlength', inf);

switch eventformat

  case 'empty'
    fprintf(1,'Pretending to write %i events...\n', length(event));

  case 'fcdc_global'
    % store it in a global variable, this is only for debugging
    if isempty(event_queue) || ~isstruct(event_queue)
      event_queue = event;
    else
      event_queue = appendstruct(event_queue, event);
    end

  case 'fcdc_rfb'
   % remote frame buffer, i.e. VNC server on another computer
   [password, host, port] = filetype_check_uri(filename);
   rfbevent(sprintf('%s:%d', host, port), password, event.type, event.value);

  case 'fcdc_buffer'
    % read from a networked buffer for realtime analysis
    [host, port] = filetype_check_uri(filename);

    if strcmp(append, 'no')
      buffer('flush_evt', [], host, port);  % flush event
    end

	% the MEX file now can handle various MATLAB types directly and respects the fields
	% sample, offset, duration
	%   -- these must all be numeric and non-empty (only first element is of interest)
	% type, value
	%   -- these can be strings or any numeric type (double, single, [u]int[8-64])
	%      will be transmitted as if vectorised

  % hack to fix empty event.offset and event.duration fields
  % Maybe should track down why they're empty? But this works for now
  % ES, 10-may-2019
  if ~isempty(event)
    for k = 1:numel(event)
      if isfield(event(k), 'offset') && isempty(event(k).offset)
        event(k).offset = 0;
      end
      if isfield(event(k), 'duration') && isempty(event(k).duration)
        event(k).duration = 0;
      end
    end
  end
	buffer('put_evt', event, host, port);

	% SK: There was some code here for firing up a FieldTrip buffer locally,
	% but this is very likely to be senseless because we have no proper header
	% information here. Please explicitly use ft_create_buffer instead.

  case 'fcdc_serial'
    % this code is moved to a separate file
    write_serial_event(filename, event);

  case 'fcdc_mysql'
    % check that the required low-level toolbox is available
    ft_hastoolbox('mysql', 1);
    % write to a MySQL server listening somewhere else on the network
    db_open(filename);
    for i=1:length(event)
      if db_blob
        % insert the structure into the database table as a binary blob
        db_insert_blob('fieldtrip.event', 'msg', event(i));
      else
        % make a structure with the same elements as the fields in the database table
        s = struct;
        % these fields also exist as elements in the table and as such can be used for filtering
        if isa(event(i).type, 'char')
          s.type = event(i).type;
        end
        if isa(event(i).value, 'numeric') && numel(event(i).value)==1
          s.value = event(i).value;
        end
        if isa(event(i).sample, 'numeric') && numel(event(i).sample)==1
          s.sample = event(i).sample;
        end
        if isa(event(i).sample, 'numeric') && numel(event(i).offset)==1
          s.offset = event(i).offset;
        end
        if isa(event(i).sample, 'numeric') && numel(event(i).duration)==1
          s.duration = event(i).duration;
        end
        % insert the structure into the database table
        db_insert('fieldtrip.event', s);
      end
    end

   case 'fcdc_fifo'

      % these are opened in blocking mode, i.e. reading/writing will block until boths sides are connected
      fifo = filetype_check_uri(filename);

      if ~exist(fifo,'file')
          ft_warning('the FIFO %s does not exist; attempting to create it', fifo);
          system(sprintf('mkfifo -m 0666 %s',fifo));
      end

      fid = fopen_or_error(fifo, 'w');
      for i=1:length(event)

        try
          % convert the event into a network message
          msg = mxSerialize(event(i));
          num = fwrite(fid, msg, 'uint8');
        catch
          ft_warning(lasterr);
        end

        if num~=length(msg)
          ft_error('problem writing to FIFO %s', fifo);
        end
      end
      fclose(fid);

    case 'fcdc_tcp'

        % TCP network socket
        [host, port] = filetype_check_uri(filename);

        con=pnet('tcpconnect',host,port);

        pnet(con,'setwritetimeout',1);

        if con~=-1,

            try % Failsafe

                for i=1:length(event)

                    % convert the event into a network message
                    msg = mxSerialize(event(i));

                    % tell the message daemon that a message will be sent, and send it
                    pnet(con,'printf',num2str(msg));
                    pnet(con,'printf','\n');
                end
%            catch
%                ft_warning(lasterr);
            end

            pnet(con,'close');
        end

    case 'fcdc_udp'

        % UDP network socket

        [host, port] = filetype_check_uri(filename);
        udp=pnet('udpsocket',port);

        if udp~=-1,
            try % Failsafe

                for i=1:length(event)

                    % convert the event into a network message
                    msg = mxSerialize(event(i));

                    % tell the message daemon that a message will be sent, and send it
                    pnet(udp,'write',uint8(msg),1000);
                    pnet(udp,'writepacket',host,port);   % Send buffer as UDP packet to host
                end

            catch
              ft_warning(lasterr);
            end
            pnet(udp,'close');
        end

    otherwise
        % assume that it is a file. Since the file probably does not yet
        % exist, determine its type by only looking at the extension
        if filetype_check_extension(filename, '.mat')
            % write the events to a MATLAB file
            if exist(filename,'file') && strcmp(append, 'yes')
                try
                    tmp = load(filename, 'event');
                    event = cat(1, tmp.event(:), event(:));
                catch
                    event = event(:);
                end
                % optionally restric the length of the event queue to flush old events
                if isfinite(maxqlength) && isreal(maxqlength) && (maxqlength>0) && (length(event)>maxqlength)
                    event = event(end-maxqlength+1:end);
                    % NOTE: this could be done using the filter event function, but
                    % then this is just a temporary solution that will probably be
                    % removed in a future versions of the code
                end
                save(filename, 'event', '-append', '-v6');
                % NOTE: the -append option in this call to the save function does
                % not actually do anything useful w.r.t. the event variable since the
                % events are being appended in the code above and the the save function
                % will just overwrite the existing event variable in the file.
                % However, if there are other variables in the file (whatever) then the
                % append option preservs them
            else
                save(filename, 'event', '-v6');
            end
        else
            ft_error('unsupported file type')
        end
end

author: spm


text:
function disp(F)
%DISP   Display a CHEBFUN2 to the command line.
%
% See also DISPLAY.

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

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

% Get display style and remove trivial empty CHEBFUN2 case.
if ( isempty(F) )
    fprintf('    empty chebfun2\n')
    if ( loose )
        fprintf('\n');
    end
    return
end

% Get information that we want to display:
dom = F.domain;                           % Domain
len = length(F);                          % Numerical rank
[xx, yy] = meshgrid(dom(1:2), dom(3:4));
vals = feval(F, xx, yy ).';               % Corner values
vals = vals(:);
vscl = vscale(F);                         % vertical scale

% Check underlying tech is a TRIGTECH:
techCol = get(F.cols.funs{1}, 'tech');
techRow = get(F.rows.funs{1}, 'tech');

% Display the information:
if ( isa(techCol(), 'trigtech') && isa(techRow(), 'trigtech') )
    disp('   chebfun2 object  (trig)')
elseif ( isa(techRow(), 'trigtech') )
    disp('   chebfun2 object  (trig in x)')
elseif ( isa(techCol(), 'trigtech') )
    disp('   chebfun2 object  (trig in y)')
else
    disp('   chebfun2 object')
end
fprintf('       domain                 rank       corner values\n');
if ( isreal(vals) )
    fprintf('[%4.2g,%4.2g] x [%4.2g,%4.2g]   %6i     [%4.2g %4.2g %4.2g %4.2g]\n', ...
        dom(1), dom(2), dom(3), dom(4), len, vals);
else
    fprintf('[%4.2g,%4.2g] x [%4.2g,%4.2g]   %6i     [  complex values  ]\n', ...
        dom(1), dom(2), dom(3) , dom(4), len);
end
fprintf('vertical scale = %3.2g \n', vscl)

if ( loose )
    fprintf('\n');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% SYNTAX:
% [pValuesS, pValuesS_cond]  = hmrS_CalcPvalue(yRuns, stimRuns, mlActRuns, baselineRange, hrfTimeWindow)
%
% UI NAME:
% Pvalues_on_Session
%
% DESCRIPTION:
% Calculate the p-value matrix for a single session.
%
% INPUTS:
% yRuns:
% stimRuns:
% mlActRuns:
% baselineRange:
% hrfTimeWindow:
%
% OUTPUTS:
% pValuesS:
% pValuesS_cond
% USAGE OPTIONS:
% Pvalues_on_Session_Concentration_Data: [pValuesS, pValuesS_cond] = hmrS_CalcPvalue(dcRuns, stimRuns, mlActRuns, baselineRange, hrfTimeWindow)
%
% PARAMETERS:
% baselineRange: [-2.0, 0.0]
% hrfTimeWindow: [-2.0, 20.0]

function [pValuesS, pValuesS_cond] = hmrS_CalcPvalue(yRuns, stimRuns, mlActRuns, baselineRange, hrfTimeWindow)

pValues = cell(length(yRuns{1}),1);

% extract fq and number of conditions from the first run:
snirf = SnirfClass(yRuns{1}, stimRuns{1});
t = snirf.GetTimeCombined();
s = snirf.GetStims(t);
ncond = size(s,2);
fq = abs(1/(t(1)-t(2)));
nDataBlks = length(yRuns{1});

%% BASELINE vs CONDITION, PAIRED T-TEST
for iBlk = 1:nDataBlks
    for cond = 1:ncond % for each condition
        for iRun = 1:length(yRuns)

            % Get stim vector by instantiating temporary SnirfClass object with this
            % function's stim argument as input, and then using the SnirfClass object's
            % GetStims method to convert stim to the s vector that this function needs.
            snirf = SnirfClass(yRuns{iRun}, stimRuns{iRun});
            t = snirf.GetTimeCombined();
            s = snirf.GetStims(t);     % stim matrix for run iRun is same for all of a run's data blocks

            % extract HRF at baselineRange and at hrfTimeWindow
            lst_stim = find(s(:,cond)==1);

            % get active measuremnt list for each run
            if isempty(mlActRuns{iRun})
                mlActRuns{iRun} = cell(length(nDataBlks),1);
            end
            ml    = yRuns{iRun}(iBlk).GetMeasListSrcDetPairs('reshape');

            if isempty(mlActRuns{iRun}{iBlk})
                mlActRuns{iRun}{iBlk} = ones(size(ml,1),1);
            end
            mlAct = mlActRuns{iRun}{iBlk}(1:size(ml,1));

            % GetDataTimeSeries() extract data in old homer2 dimensions (time X Hb X channel)
            y = yRuns{iRun}(iBlk).GetDataTimeSeries('reshape'); % yRuns{iRun}(iBlk).GetDataTimeSeries('reshape');  % data matrix for run iRun, data block iBlk
            for hb = 1:3 % across HbO/HbR/HbT
                for iTrial = 1:size(lst_stim,1) % across trials
                    % Hb: # of time points X # of trials X # of runs X # of channels X HbO/HbR
                    Hb_baseline(:,iTrial,iRun,:,hb,cond) = squeeze(y([lst_stim(iTrial) - round(abs(baselineRange(1))*fq)]:lst_stim(iTrial),hb,:)); % get each trial
                    Hb_peak(:,iTrial,iRun,:,hb,cond) = squeeze(y([lst_stim(iTrial) + round(abs(hrfTimeWindow(1))*fq)]:[lst_stim(iTrial) + round(hrfTimeWindow(2)*fq)],hb,:)); % get each trial

                end
            end
        end
    end

    % put together trials from all runs:
    Hb_baseline_rs = reshape(Hb_baseline, size(Hb_baseline,1),  size(Hb_baseline,2)* size(Hb_baseline,3),  size(Hb_baseline,4),size(Hb_baseline,5),size(Hb_peak,6));
    Hb_peak_rs = reshape(Hb_peak, size(Hb_peak,1),  size(Hb_peak,2)* size(Hb_peak,3),  size(Hb_peak,4),size(Hb_peak,5),size(Hb_peak,6));

    % take the mean in time ranges: baselineRange, hrfTimeWindow
    for hb = 1:3 % HbO/HbR/HbT
        for ch=1:size(Hb_peak_rs,3) % across channels
            MEAN_Hb_baseline(:,ch,hb)= nanmean(squeeze(Hb_baseline_rs(:,:,ch,hb)),1);
            MEAN_Hb_peak(:,ch,hb)= nanmean(squeeze(Hb_peak_rs(:,:,ch,hb)),1);
        end
    end

    % get stats
    if isempty(pValues{iBlk})
        for ch = 1:size(MEAN_Hb_peak,2) % channels
                 for hb = 1:3% HbO/HbR/HbT
                    [h,p,c,stats] = ttest(MEAN_Hb_baseline(:,ch,hb),(MEAN_Hb_peak(:,ch,hb)));
                    pValuesS(hb,ch,cond) = p;

                    % 2) Pvalue SHOULD BE CONVERTED TO DATA CLASS IN THE RIGHT DIM AT THE END

                end
        end
    end
end

pValuesS(:,find(mlAct==0),:) = NaN;

%% CONDITION vs CONDITION, UNPAIRED T-TEST
% get all combinations of conditions
cond_2_comb = sort(combnk(1:ncond,2));
% extract HRF at hrfTimeWindow from the cond combination
for i = 1:ncond
    lst_stim_all{i} = find(s(:,i)==1);
end

for iBlk = 1:length(nDataBlks)
    for comb_inx = 1:size(cond_2_comb,1) % for each condition
        % get current combin.
        foo = cond_2_comb(comb_inx,:);

        for iRun = 1:length(yRuns)
            % get active measuremnt list for each run
            if isempty(mlActRuns{iRun})
                mlActRuns{iRun} = cell(length(nDataBlks),1);
            end
            ml = yRuns{iRun}(iBlk).GetMeasListSrcDetPairs('reshape');

            if isempty(mlActRuns{iRun}{iBlk})
                mlActRuns{iRun}{iBlk} = ones(size(ml,1),1);
            end
            mlAct = mlActRuns{iRun}{iBlk}(1:size(ml,1));

            % GetDataTimeSeries() extract data in old homer2 dimensions (time X Hb X channel)
            y = yRuns{iRun}(iBlk).GetDataTimeSeries('reshape');  % data matrix for run iRun, data block iBlk
            for hb = 1:3 % across HbO/HbR/HbT
                for iTrial = 1:size(lst_stim_all{foo(1)},1) % across trials
                    % Hb: # of time points X # of trials X # of runs X # of channels X HbO/HbR
                    Hb_peak1(:,iTrial,iRun,:,hb) = squeeze(y([lst_stim_all{foo(1)}(iTrial) + round(abs(hrfTimeWindow(1))*fq)]:[lst_stim_all{foo(1)}(iTrial) + round(hrfTimeWindow(2)*fq)],hb,:)); % get each trial
                end
                for iTrial = 1:size(lst_stim_all{foo(2)},1) % across trials
                    % Hb: # of time points X # of trials X # of runs X # of channels X HbO/HbR
                    Hb_peak2(:,iTrial,iRun,:,hb) = squeeze(y([lst_stim_all{foo(2)}(iTrial) + round(abs(hrfTimeWindow(1))*fq)]:[lst_stim_all{foo(2)}(iTrial) + round(hrfTimeWindow(2)*fq)],hb,:)); % get each trial
                end
            end
        end
    end

    % put together trials from all runs:
    Hb_peak_rs1 = reshape(Hb_peak1, size(Hb_peak1,1),  size(Hb_peak1,2)* size(Hb_peak1,3),  size(Hb_peak1,4),size(Hb_peak1,5));
    Hb_peak_rs2 = reshape(Hb_peak2, size(Hb_peak2,1),  size(Hb_peak2,2)* size(Hb_peak2,3),  size(Hb_peak2,4),size(Hb_peak2,5));

    % take the mean in time range hrfTimeWindow
    for hb = 1:3 % HbO/HbR/HbT
        for ch=1:size(Hb_peak_rs,3) % across channels
            MEAN_Hb_peak1(:,ch,hb)= nanmean(squeeze(Hb_peak_rs1(:,:,ch,hb)),1);
            MEAN_Hb_peak2(:,ch,hb)= nanmean(squeeze(Hb_peak_rs2(:,:,ch,hb)),1);

        end
    end

    % get stats
    if ~exist('pValuesS_cond')
        for ch = 1:size(MEAN_Hb_peak,2) % channels
                for hb = 1:3% HbO/HbR/HbT
                    [h,p,c,stats] = ttest2(MEAN_Hb_peak1(:,ch,hb),(MEAN_Hb_peak2(:,ch,hb)));
                    pValuesS_cond(hb,ch,comb_inx) = p;
                    % or                     pValuesS_cond(foo(1),foo(2),hb,ch) = p;

                    % 2) Pvalue SHOULD BE CONVERTED TO DATA CLASS IN THE RIGHT DIM AT THE END

                end
        end
    end
end
pValuesS_cond(:,find(mlAct==0),:) = NaN;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: BUNPC


text:
function flag = areBinaryWeights(weights)
% Check the input vector is (approximately) a binary vectory
flag = true;
th = 1e-12;
for i = length(weights):-1:1
    if weights(i) > th && abs(1 - weights(i)) > th
        flag = false;
        break
    end
end
end

author: MIT-SPARK

repo: GNC-and-ADAPT


text:
function i4vec_min_index_test ( )

%*****************************************************************************80
%
%% I4VEC_MIN_INDEX_TEST tests I4VEC_MIN_INDEX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_MIN_INDEX_TEST\n' );
  fprintf ( 1, '  For an integer vector:\n' );
  fprintf ( 1, '  I4VEC_MIN_INDEX:          a minimal index;\n' );

  seed = 123456789;
  b = -n;
  c = n;

  [ a, seed ] = i4vec_uniform_ab ( n, b, c, seed );

  i4vec_print ( n, a, '  Input vector:' );

  fprintf ( 1, '\n' );

  ival = i4vec_min_index ( n, a );
  fprintf ( 1, '  Minimum index: %d\n', ival );

  return
end


text:
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %   PULSED SPOTLIGHT SAR SIMULATION AND RECONSTRUCTION   %
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

colormap(gray(256))
cj=sqrt(-1);
pi2=2*pi;
%
c=3e8;                   % propagation speed
f0=50e6;                 % baseband bandwidth is 2*f0
w0=pi2*f0;
fc=200e6;                % carrier frequency
wc=pi2*fc;
lambda_min=c/(fc+f0);    % Wavelength at highest frequency
lambda_max=c/(fc-f0);    % Wavelength at lowest frequency
kc=(pi2*fc)/c;           % wavenumber at carrier frequency
kmin=(pi2*(fc-f0))/c;    % wavenumber at lowest frequency
kmax=(pi2*(fc+f0))/c;    % wavenumber at highest frequency
%
Xc=1000;                 % Range distance to center of target area
X0=20;                   % target area in range is within [Xc-X0,Xc+X0]
Yc=300;                  % Cross-range distance to center of target area
Y0=60;                  % target area in cross-range is within
                         % [Yc-Y0,Yc+Y0]

% Case 1: L < Y0; requires zero-padding of SAR signal in synthetic
% aperture domain
%
  L=100;                 % synthetic aperture is 2*L

% Case 2: L > Y0; slow-time Doppler subsampling of SAR signal spectrum
% reduces computation
%
% L=400;                 % synthetic aperture is 2*L

theta_c=atan(Yc/Xc);     % Squint angle
Rc=sqrt(Xc^2+Yc^2);      % Squint radial range
L_min=max(Y0,L);         % Zero-padded aperture is 2*L_min

%
Xcc=Xc/(cos(theta_c)^2); % redefine Xc by Xcc for squint processing

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% u domain parameters and arrays for compressed SAR signal %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
duc=(Xcc*lambda_min)/(4*Y0);      % sample spacing in aperture domain
                                  % for compressed SAR signal
duc=duc/1.2;                      % 10 percent guard band; this guard band
                                  % would not be sufficient for targets
                                  % outside digital spotlight filter (use
                                  % a larger guard band, i.e., PRF)
mc=2*ceil(L_min/duc);             % number of samples on aperture
uc=duc*(-mc/2:mc/2-1);            % synthetic aperture array
dkuc=pi2/(mc*duc);                % sample spacing in ku domain
kuc=dkuc*(-mc/2:mc/2-1);          % kuc array
%
dku=dkuc;                         % sample spacing in ku domain

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%    u domain parameters and arrays for SAR signal     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if Yc-Y0-L < 0,                            % minimum aspect angle
 theta_min=atan((Yc-Y0-L)/(Xc-X0));
else,
 theta_min=atan((Yc-Y0-L)/(Xc+X0));
end;
theta_max=atan((Yc+Y0+L)/(Xc-X0));         % maximum aspect angle
%
du=pi/(kmax*(sin(theta_max)- ...
                     sin(theta_min))); % sample spacing in aperture
                                       % domain for SAR signal
du=du/1.4;                        % 20 percent guard band
m=2*ceil(pi/(du*dku));            % number of samples on aperture
du=pi2/(m*dku);                   % readjust du
u=du*(-m/2:m/2-1);                % synthetic aperture array
ku=dku*(-m/2:m/2-1);              % ku array

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%       Fast-time domain parmeters and arrays          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
Tp=2.5e-7;                     % Chirp pulse duration
alpha=w0/Tp;                   % Chirp rate
wcm=wc-alpha*Tp;               % Modified chirp carrier
%
if Yc-Y0-L < 0,
 Rmin=Xc-X0;
else,
 Rmin=sqrt((Xc-X0)^2+(Yc-Y0-L)^2);
end;
Ts=(2/c)*Rmin;                 % start time of sampling
Rmax=sqrt((Xc+X0)^2+(Yc+Y0+L)^2);
Tf=(2/c)*Rmax+Tp;              % end time of sampling
T=Tf-Ts;                       % fast-time interval of measurement
Ts=Ts-.1*T;                    % start slightly earlier (10% guard band)
Tf=Tf+.1*T;                    % end slightly later (10% guard band)
T=Tf-Ts;
Tmin=max(T,(4*X0)/(c*cos(theta_max)));  % Minimum required T
%
dt=1/(4*f0);                 % Time domain sampling (guard band factor 2)
n=2*ceil((.5*Tmin)/dt);      % number of time samples
t=Ts+(0:n-1)*dt;             % time array for data acquisition
dw=pi2/(n*dt);               % Frequency domain sampling
w=wc+dw*(-n/2:n/2-1);        % Frequency array (centered at carrier)
k=w/c;                       % Wavenumber array
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Resolution for Broadside: (x,y) domain rotated by theta_c %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

DX=c/(4*f0);                      % range resolution (broadside)
DY=(Xcc*lambda_max)/(4*L);         % cross-range resolution (broadside)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%           Parameters of Targets                 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
ntarget=9;                        % number of targets
% Set ntarget=1 to see "clean" PSF of target at origin
% Try this with other targets

% xn: range;            yn= cross-range;    fn: reflectivity
  xn=zeros(1,ntarget);  yn=xn;              fn=xn;

% Targets within digital spotlight filter
%
  xn(1)=0;              yn(1)=0;            fn(1)=1;
  xn(2)=.7*X0;          yn(2)=-.6*Y0;       fn(2)=1.4;
  xn(3)=0;              yn(3)=-.85*Y0;      fn(3)=.8;
  xn(4)=-.5*X0;         yn(4)=.75*Y0;       fn(4)=1.;
  xn(5)=-.5*X0+DX;      yn(5)=.75*Y0+DY;    fn(5)=1.;

% Targets outside digital spotlight filter
% (Run the code with and without these targets)
%
  xn(6)=-1.2*X0;        yn(6)=.75*Y0;       fn(6)=1.;
  xn(7)=.5*X0;          yn(7)=1.25*Y0;      fn(7)=1.;
  xn(8)=1.1*X0;         yn(8)=-1.1*Y0;      fn(8)=1.;
  xn(9)=-1.2*X0;        yn(9)=-1.75*Y0;     fn(9)=1.;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                   SIMULATION                    %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
s=zeros(n,mc);     % SAR signal array
%
for i=1:ntarget;   % Loop for each target
 td=t(:)*ones(1,mc)-2*ones(n,1)*sqrt((Xc+xn(i)).^2+(Yc+yn(i)-uc).^2)/c;
 s=s+fn(i)*exp(cj*wcm*td+cj*alpha*(td.^2)).*(td >= 0 & td <= Tp & ...
   ones(n,1)*abs(uc) <= L & t(:)*ones(1,mc) < Tf);
end;
%
s=s.*exp(-cj*wc*t(:)*ones(1,mc));      % Fast-time baseband conversion

% User may apply a slow-time domain window, e.g., power window, on
% simulated SAR signal array "s" here.

G=abs(s)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(t,uc,256-cg*(G-ng));
axis('square');axis('xy')
xlabel('Fast-time t, sec')
ylabel('Synthetic Aperture (Slow-time) U, meters')
title('Measured Spotlight SAR Signal')
print P5.1.ps
pause(1)
%

td0=t(:)-2*sqrt(Xc^2+Yc^2)/c;
s0=exp(cj*wcm*td0+cj*alpha*(td0.^2)).*(td0 >= 0 & td0 <= Tp);
s0=s0.*exp(-cj*wc*t(:));            % Baseband reference fast-time signal

s=ftx(s).*(conj(ftx(s0))*ones(1,mc));  % Fast-time matched filtering
%
G=abs(iftx(s))';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
tm=(2*Rc/c)+dt*(-n/2:n/2-1);    % fast-time array after matched filtering
image(tm,uc,256-cg*(G-ng));
axis('square');axis('xy')
xlabel('Fast-time t, sec')
ylabel('Synthetic Aperture (Slow-time) U, meters')
title('SAR Signal after Fast-time Matched Filtering')
print P5.2.ps
pause(1)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Slow-time baseband conversion for squint %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
kus=2*kc*sin(theta_c)*ones(1,n);     % Doppler frequency shift in ku
                                     % domain due to squint
%
s=s.*exp(-cj*kus(:)*uc);             % slow-time baseband conversion
fs=fty(s);

% Display aliased SAR spectrum
%
G=abs(fs)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(k*c/pi2,kuc,256-cg*(G-ng));
axis('square');axis('xy')
xlabel('Fast-time Frequency, Hertz')
ylabel('Synthetic Aperture (Slow-time) Frequency Ku, rad/m')
title('Aliased Spotlight SAR Signal Spectrum')
print P5.3.ps
pause(1)

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Digital Spotlighting and Bandwidth Expansion in ku Domain  %%
%%          via Slow-time Compression and Decompression        %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
s=s.*exp(cj*kus(:)*uc);      % Original signal before baseband
                             % conversion for squint

cs=s.*exp(cj*2*(k(:)*ones(1,mc)).* ...
 (ones(n,1)*sqrt(Xc^2+(Yc-uc).^2))-cj*2*k(:)*Rc*ones(1,mc));% compression
fcs=fty(cs);            % F.T. of compressed signal w.r.t. u
%
G=abs(fcs)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(k*c/pi2,kuc,256-cg*(G-ng));
axis('square');axis('xy')
xlabel('Fast-time Frequency, Hertz')
ylabel('Synthetic Aperture (Slow-time) Frequency Ku, rad/m')
title('Compressed Spotlight SAR Signal Spectrum')
print P5.4.ps
pause(1)
%
fp=iftx(fty(cs));      % Narrow-bandwidth Polar Format Processed
                       % reconstruction
%
PH=asin(kuc/(2*kc));   % angular Doppler domain
R=(c*tm)/2;            % range domain mapped from reference
                       % fast-time domain
%
% Full Aperture Digital-Spotlight Filter
%
W_d=((abs(R(:)*cos(PH+theta_c)-Xc) < X0).* ...
    (abs(R(:)*sin(PH+theta_c)-Yc) < Y0));
%
G=(abs(fp)/max(max(abs(fp)))+.1*W_d)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image((Rc/Xc)*(.5*c*tm-Rc),(kuc*Rc)/(2*kc),256-cg*(G-ng));
xlabel('Range x, m')
ylabel('Cross-range y, m')
title('Polar Format SAR Reconstruction with Digital Spotlight Filter')
axis image; axis xy;
print P5.5.ps
pause(1)

fd=fp.*W_d;                % Digital Spotlight Filtering
fcs=ftx(fd);               % Transform to (omega,ku) domain

% Zero-padding in ku domain for slow-time upsampling
%
mz=m-mc;        % number is zeros
fcs=(m/mc)*[zeros(n,mz/2),fcs,zeros(n,mz/2)];
%
cs=ifty(fcs);              % Transform to (omega,u) domain

s=cs.*exp(-cj*2*(k(:)*ones(1,m)).* ...
 (ones(n,1)*sqrt(Xc^2+(Yc-u).^2))+cj*2*k(:)*Rc*ones(1,m));% decompression

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           CAUTION                             %
% For TDC or backprojection, do not subsample in Doppler domain %
% and do not perform slow-time baseband conversion               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
s_ds=s;                    % Save s(omega,u) array for TDC and
                           % backprojection algorithms

%
s=s.*exp(-cj*kus(:)*u);    % Slow-time baseband conversion for squint
fs=fty(s);                 % Digitally-spotlighted SAR signal spectrum
%
G=abs(fs)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(k*c/pi2,ku,256-cg*(G-ng));
axis('square');axis('xy')
xlabel('Fast-time Frequency, Hertz')
ylabel('Synthetic Aperture (Slow-time) Frequency Ku, rad/m')
title('Spotlight SAR Signal Spectrum after DS & Upsampling')
print P5.6.ps
pause(1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%    SLOW-TIME DOPPLER SUBSAMPLING     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if Y0 < L,
 ny=2*ceil(1.2*Y0/du);      % Number of samples in y domain
                            % 20 percent guard band
 ms=floor(m/ny);            % subsampling ratio
 tt=floor(m/(2*ms));
 I=m/2+1-tt*ms:ms:m/2+1+(tt-1)*ms; % subsampled index in ku domain
 [tt,ny]=size(I);           % number of subsamples
 fs=fs(:,I);                % subsampled SAR signal spectrum
 ky=ku(I);                  % subsampled ky array
 dky=dku*ms;                % ky domain sample spacing
else,
 dky=dku;
 ny=m;
 ky=ku;
end;

dy=pi2/(ny*dky);            % y domain sample spacing
y=dy*(-ny/2:ny/2-1);        % cross-range array

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%             RECONSTRUCTION           %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
ky=ones(n,1)*ky+kus(:)*ones(1,ny);       % ky array
kx=(4*k(:).^2)*ones(1,ny)-ky.^2;
kx=sqrt(kx.*(kx > 0));                  % kx array
%
plot(kx(1:20:n*ny),ky(1:20:n*ny),'.')
xlabel('Spatial Frequency k_x, rad/m')
ylabel('Spatial Frequency k_y, rad/m')
title('Spotlight SAR Spatial Frequency Data Coverage')
axis image; axis xy
print P5.7.ps
pause(1)
%
kxmin=min(min(kx));
kxmax=max(max(kx));
dkx=pi/X0;        % Nyquist sample spacing in kx domain
nx=2*ceil((.5*(kxmax-kxmin))/dkx); % Required number of
                      % samples in kx domain;
                      % This value will be increased slightly
                      % to avoid negative array index
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                                                         %%%
%%%   FIRST TWO OPTIONS FOR RECONSTRUCTION:                 %%%
%%%                                                         %%%
%%%     1. 2D Fourier Matched Filtering and Interpolation   %%%
%%%     2. Range Stacking                                   %%%
%%%                                                         %%%
%%%     Note: For "Range Stacking," make sure that the      %%%
%%%           arrays nx, x, and kx are defined.             %%%
%%%                                                         %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%    2D FOURIER MATCHED FILTERING AND INTERPOLATION    %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Matched Filtering
%
fs0=(kx > 0).*exp(cj*kx*Xc+cj*ky*Yc+cj*.25*pi ...
           -cj*2*k(:)*ones(1,ny)*Rc); % reference signal complex conjugate
fsm=fs.*fs0;     % 2D Matched filtering

% Interpolation
%
is=8;       % number of neighbors (sidelobes) used for sinc interpolator
I=2*is+1;
kxs=is*dkx; % plus/minus size of interpolation neighborhood in KX domain
%
nx=nx+2*is+4;  % increase number of samples to avoid negative
               %  array index during interpolation in kx domain
KX=kxmin+(-is-2:nx-is-3)*dkx;     % uniformly-spaced kx points where
                                  % interpolation is done
kxc=KX(nx/2+1);                   % carrier frequency in kx domain
KX=KX(:)*ones(1,ny);
%
F=zeros(nx,ny);         % initialize F(kx,ky) array for interpolation

for i=1:n;                       % for each k loop
  i                              % print i to show that it is running
 icKX=round((kx(i,:)-KX(1,1))/dkx)+1; % closest grid point in KX domain
 cKX=KX(1,1)+(icKX-1)*dkx;            % and its KX value
 ikx=ones(I,1)*icKX+[-is:is]'*ones(1,ny);
 ikx=ikx+nx*ones(I,1)*[0:ny-1];
 nKX=KX(ikx);
 SINC=sinc((nKX-ones(I,1)*kx(i,:))/dkx);             % interpolating sinc
 HAM=.54+.46*cos((pi/kxs)*(nKX-ones(I,1)*kx(i,:)));  % Hamming window
         %%%%%   Sinc Convolution (interpolation) follows  %%%%%%%%
 F(ikx)=F(ikx)+(ones(I,1)*fsm(i,:)).*(SINC.*HAM);
end
%
%  DISPLAY interpolated spatial frequency domain image F(kx,ky)

KX=KX(:,1).';
KY=ky(1,:);

G=abs(F)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(KX,KY+kus(1),256-cg*(G-ng));
axis image; axis xy
xlabel('Spatial Frequency k_x, rad/m')
ylabel('Spatial Frequency k_y, rad/m')
title('Wavefront Spotlight SAR Reconstruction Spectrum')
print P5.8.ps
pause(1)

%
f=iftx(ifty(F));     % Inverse 2D FFT for spatial domain image f(x,y)
%
dx=pi2/(nx*dkx);     % range sample spacing in reconstructed image
x=dx*(-nx/2:nx/2-1); % range array
%
% Display SAR reconstructed image

G=abs(f)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(Xc+x,Yc+y,256-cg*(G-ng));axis([Xc-X0 Xc+X0 Yc-Y0 Yc+Y0]);
axis image; axis xy
xlabel('Range X, meters')
ylabel('Cross-range Y, meters')
title('Wavefront Spotlight SAR Reconstruction')
print P5.9.ps
pause(1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%   SAR Image Compression (for Spotlight System)  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
Fc=ftx(fty(f.* ...
  exp(cj*kxc*x(:)*ones(1,ny)+cj*ones(nx,1)*2*kc*sin(theta_c)*y ...
 -cj*2*kc*sqrt(((Xc+x(:)).^2)*ones(1,ny)+ones(nx,1)*((Yc+y).^2)))));
G=abs(Fc)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(KX,KY+kus(1),256-cg*(G-ng));
axis image; axis xy
xlabel('Spatial Frequency k_x, rad/m')
ylabel('Spatial Frequency k_y, rad/m')
title('Compressed Spotlight SAR Reconstruction Spectrum')
print P5.10.ps
pause(1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%      RANGE STACK WAVEFRONT RECONSTRUCTION       %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
f_stack=zeros(nx,ny); % Initialize reconstruction array in (x,y) domain
for i=1:nx; i        % Stack's loop for reconstruction at each range
 f_stack(i,:)=ifty(sum(fs.*exp(cj*kx*(Xc+x(i))+cj*ky*Yc ...
  +cj*.25*pi-cj*2*k(:)*ones(1,ny)*Rc)));
end;

% Remove carrier in range domain
f_stack=f_stack.*exp(-cj*x(:)*kxc*ones(1,ny));
%
f_stack=f_stack/nx; % Scale it for comparison with Fourier interpolation
                    % Use "f_stack-f" to display difference of two
                    % reconstructions

G=abs(f_stack)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(Xc+x,Yc+y,256-cg*(G-ng));axis([Xc-X0 Xc+X0 Yc-Y0 Yc+Y0]);
axis image; axis xy
xlabel('Range X, meters')
ylabel('Cross-range Y, meters')
title('Range Stack Spotlight SAR Reconstruction')
print P5.11.ps
pause(1)

F_stack=ftx(fty(f_stack)); % Reconstruction array in spatial frequency
                           % domain; Use "F_stack-F" to display
                           %  difference of two reconstructions
%
G=abs(F_stack)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(KX,KY+kus(1),256-cg*(G-ng));
axis image; axis xy
xlabel('Spatial Frequency k_x, rad/m')
ylabel('Spatial Frequency k_y, rad/m')
title('Range Stack Spotlight SAR Reconstruction Spectrum')
print P5.12.ps
pause(1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%     TIME DOMAIN CORRELATION RECONSTRUCTION      %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
f_tdc=zeros(nx,ny); % Initialize reconstruction array in (x,y) domain

for i=1:nx; i
   for j=1:ny;
      t_ij=(2*sqrt((x(i)+Xc)^2+(y(j)+Yc-u).^2))/c;
      f_tdc(i,j)=sum(sum(s_ds.*exp(cj*w(:)*(t_ij-tm(n/2+1))).* ...
         (ones(n,1)*(t_ij >= Ts & t_ij <= Tf))));
   end;
end;

% Remove carrier in range domain
f_tdc=f_tdc.*exp(-cj*x(:)*kxc*ones(1,ny));
%

% Remove carrier in cross-range domain (squint mode)
f_tdc=f_tdc.*exp(-cj*ones(nx,1)*2*kc*sin(theta_c)*y);

G=abs(f_tdc)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(Xc+x,Yc+y,256-cg*(G-ng));axis([Xc-X0 Xc+X0 Yc-Y0 Yc+Y0]);
axis image; axis xy
xlabel('Range X, meters')
ylabel('Cross-range Y, meters')
title('TDC Spotlight SAR Reconstruction')
print P5.13.ps
pause(1)

F_tdc=ftx(fty(f_tdc));     % Reconstruction array in spatial frequency
%
G=abs(F_tdc)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(KX,KY+kus(1),256-cg*(G-ng));
axis image; axis xy
xlabel('Spatial Frequency k_x, rad/m')
ylabel('Spatial Frequency k_y, rad/m')
title('TDC Spotlight SAR Reconstruction Spectrum')
print P5.14.ps
pause(1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%           BACKPROJECTION RECONSTRUCTION         %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
f_back=zeros(nx,ny); % Initialize reconstruction array in (x,y) domain
n_ratio=100;              % Upsampling ratio in fast-time domain
nu=n_ratio*n;             % Size of upsampled s(t,u) array in t domain
nz=nu-n;                  % Number of zeros
dtu=(n/nu)*dt;            % Fast-time sample spacing of upsampled array
tu=dtu*(-nu/2:nu/2-1);    % Upsampled reference fast-time array
X=x(:)*ones(1,ny);
Y=ones(nx,1)*y;

for j=1:m; j
   t_ij=(2*sqrt((X+Xc).^2+(Y+Yc-u(j)).^2))/c;
   t_ij=round((t_ij-tm(n/2+1))/dtu)+nu/2+1;
   it_ij=(t_ij > 0 & t_ij <= nu);
   t_ij=t_ij.*it_ij+nu*(1-it_ij);
   S=ifty([zeros(1,nz/2),s_ds(:,j).',zeros(1,nz/2)])...
      .*exp(cj*wc*tu);
   S(nu)=0;
   f_back=f_back+S(t_ij);
end;

clear X Y

% Remove carrier in range domain
f_back=f_back.*exp(-cj*x(:)*kxc*ones(1,ny));
%

% Remove carrier in cross-range domain (squint mode)
f_back=f_back.*exp(-cj*ones(nx,1)*2*kc*sin(theta_c)*y);

G=abs(f_back)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(Xc+x,Yc+y,256-cg*(G-ng));axis([Xc-X0 Xc+X0 Yc-Y0 Yc+Y0]);
axis image; axis xy
xlabel('Range X, meters')
ylabel('Cross-range Y, meters')
title('Backprojection Spotlight SAR Reconstruction')
print P5.15.ps
pause(1)

F_back=ftx(fty(f_back));     % Reconstruction array in spatial frequency
%
G=abs(F_back)';
xg=max(max(G)); ng=min(min(G)); cg=255/(xg-ng);
image(KX,KY+kus(1),256-cg*(G-ng));
axis image; axis xy
xlabel('Spatial Frequency k_x, rad/m')
ylabel('Spatial Frequency k_y, rad/m')
title('Backprojection Spotlight SAR Reconstruction Spectrum')
print P5.16.ps
pause(1)


text:
function [ alphas, betas, scaling, finalLikelihood] = CCRF_training_gradient_descent( nIterations, nExamples, learningRate, threshold, x, y, yUnnormed, masks, alphas, betas, lambda_a, lambda_b, similarityFNs, useIndicators, verbose)
%GRADIENTDESCENTCCRF Performs CCRF gradient descen given the initial state
%and gradient descent parameters
%   Detailed explanation goes here

    if(verbose)
        logLikelihood = zeros(round(nIterations/10)+1, 1);
        alphaTrack = zeros(nIterations, numel(alphas));
        betaTrack = zeros(nIterations, numel(betas));
    end

    logAlphas = log(alphas);
    logBetas = log(betas);

    K = numel(similarityFNs);

    %calculate similarity measures for each of the sequences
    Similarities = cell(nExamples, 1);
    PrecalcQ2s = cell(nExamples,1);
    PrecalcQ2sFlat = cell(nExamples,1);

    PrecalcYqDs = zeros(nExamples, K);

    for q = 1 : nExamples

        yq = y{q};
        xq = x{q};
        mask = masks{q};

        n = size(yq, 1);
        Similarities{q} = zeros([n, n, K]);
%         PrecalcQ2s{q} = zeros([n, n, K]);
        PrecalcQ2s{q} = cell(K,1);
%         PrecalcQ2sFlat{q} = cell(K,1);
        PrecalcQ2sFlat{q} = zeros((n*(n+1))/2,K);
        % go over all of the similarity metrics and construct the
        % similarity matrices
        for k=1:K
            Similarities{q}(:,:,k) = similarityFNs{k}(xq, mask);
            S = Similarities{q}(:,:,k);
            D =  diag(sum(S));
            B = D - S;
%             PrecalcQ2s{q}(:,:,k) = B;
            PrecalcQ2s{q}{k} = B;
%             PrecalcQ2sFlat{q}{k} = PrecalcQ2s{q}{k}(logical(tril(ones(size(S)))));
            PrecalcQ2sFlat{q}(:,k) = B(logical(tril(ones(size(S)))));
            PrecalcYqDs(q,k) = -yq'*B*yq;
        end
    end

    %stochastic gradient descent
    for iter = 1 : nIterations
        prevAlphas = alphas;
        prevBetas = betas;

        for q = 1 : nExamples

            yq = y{q};
            xq = x{q};
            mask = masks{q};

            PrecalcQ2 = PrecalcQ2s{q};
            PrecalcQ2Flat = PrecalcQ2sFlat{q};
            [ logGradientsAlphas, logGradientsBetas] = gradientCCRF(alphas, betas, lambda_a, lambda_b, PrecalcQ2, xq, yq, mask, PrecalcYqDs(q, :), useIndicators, PrecalcQ2Flat);

%             [logGradientAlphasAnalytical, logGradientBetasAnalytical] = gradientAnalytical(PrecalcQ2, alphas, betas, lambda, xq, yq, mask);
%
%             diffInGradientsAlpha = mean(abs(logGradientsAlphas - logGradientAlphasAnalytical));
%             diffInGradientsBeta = mean(abs(logGradientsBetas - logGradientBetasAnalytical));

            %update log alpha
            logAlphas = logAlphas + learningRate * logGradientsAlphas;
            alphas = exp(logAlphas);

            %update log beta
            logBetas = logBetas + learningRate * logGradientsBetas;
            betas = exp(logBetas);

            if(verbose)
                %record alpha and beta values for each iteration for debug purposes
                alphaTrack(iter,:) = alphas(:);
                betaTrack(iter,:) = betas;
            end
        end

        %check for convergence
        if (norm([prevAlphas;prevBetas] - [alphas;betas])/norm([prevAlphas;prevBetas]) < threshold || norm([logGradientsAlphas;logGradientsBetas]) < threshold)
            break;
        end

        if(verbose)
            if(mod(iter, 10)==0)
                logLikelihood(iter/10 + 1) = LogLikelihoodCCRF(y, x, masks, alphas, betas, lambda_a, lambda_b, PrecalcQ2sFlat, useIndicators);
                fprintf('Iteration %d; logL %f\n', iter, logLikelihood(iter/10 + 1));
            end

        end
    end

    % establish the scaling
    scaling = getScaling(alphas, betas, x, yUnnormed, masks, PrecalcQ2s, useIndicators);

    if(verbose)
        figure
        subplot(1,3,1)
        plot(betaTrack(1:iter,:));
        title('beta');
        subplot(1,3,2)
        plot(alphaTrack(1:iter,:))
        title('alpha');
        subplot(1,3,3)
        plot(logLikelihood(1:round(iter/10),:))
        title('log likelihood');
        finalLikelihood = LogLikelihoodCCRF(y, x, masks, alphas, betas, lambda_a, lambda_b, PrecalcQ2sFlat, useIndicators);
        fprintf('Final log likelihood at iteration %d; logL %f, learning rate %f\n', iter, finalLikelihood, learningRate);
    else
        finalLikelihood = LogLikelihoodCCRF(y, x, masks, alphas, betas, lambda_a, lambda_b, PrecalcQ2sFlat, useIndicators);
        fprintf('Final log likelihood at iteration %d; logL %f, learning rate %f\n', iter, finalLikelihood, learningRate);
    end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function [revisedModel,gapfilledReactions,replacedReactions]=debugModel(model,testResults,inputDataFolder,infoFilePath,microbeID,biomassReaction)
% This function runs a suite of debugging functions on a refined
% reconstruction produced by the DEMETER pipeline. Tests
% are performed whether or not the models can produce biomass aerobically
% and anaerobically, and whether or not unrealistically high ATP is
% produced on a complex medium.
%
% USAGE:
%
%   [revisedModel,gapfilledReactions,replacedReactions]=debugModel(model,testResults,inputDataFolder,infoFilePath,microbeID,biomassReaction)
%
% INPUTS
% model:                 COBRA model structure
% testResults:           Structure with results of test run
% inputDataFolder:       Folder with input tables with experimental data
%                        and databases that inform the refinement process
% infoFilePath           File with information on reconstructions to refine
% microbeID:             ID of the reconstructed microbe that serves as
%                        the reconstruction name and to identify it in
%                        input tables
% biomassReaction:       Reaction ID of the biomass objective function
%
% OUTPUT
% revisedModel:          Gapfilled COBRA model structure
% gapfilledReactions:    Reactions gapfilled to enable flux
% replacedReactions:     Reactions replaced because they were causing
%                        futile cycles
%
% .. Author:
%       - Almut Heinken, 09/2020

gapfilledReactions = {};
cntGF=1;
replacedReactions = {};

tol=0.0000001;

model=changeObjective(model,biomassReaction);

% implement complex medium
constraints = readtable('ComplexMedium.txt', 'Delimiter', 'tab');
constraints=table2cell(constraints);
constraints=cellstr(string(constraints));

% Load reaction and metabolite database
database=loadVMHDatabase;

% create a temporary summary file of the performed refinement
summary=struct;
summary.condGF={};
summary.targetGF={};
summary.relaxGF={};

% rebuild model
model = rebuildModel(model,database,biomassReaction);

[AerobicGrowth, AnaerobicGrowth] = testGrowth(model, biomassReaction);
if AnaerobicGrowth(1,1) < tol
    % find reactions that are preventing the model from growing
    % anaerobically
    % first gapfilling specialized for anaerobic growth
    [model,oxGapfillRxns,anaerGrowthOK] = anaerobicGrowthGapfill(model, biomassReaction, database);
    if ~isempty(oxGapfillRxns)
        summary.condGF=union(summary.condGF,oxGapfillRxns);

        gapfilledReactions{cntGF,1}=microbeID;
        gapfilledReactions{cntGF,2}='Enabling anaerobic growth';
        gapfilledReactions{cntGF,3}='Condition-specific gapfilling';
        gapfilledReactions(cntGF,4:length(oxGapfillRxns)+3)=oxGapfillRxns;
        cntGF=cntGF+1;
    end
    % then less targeted gapfilling
    model=changeRxnBounds(model,'EX_o2(e)',0,'l');
    [model,condGF,targetGF,relaxGF] = runGapfillingFunctions(model,biomassReaction,biomassReaction,'max',database,1);
    model=changeRxnBounds(model,'EX_o2(e)',-10,'l');
    % export the gapfilled reactions
    if ~isempty(condGF)
        summary.condGF=union(summary.condGF,condGF);

        gapfilledReactions{cntGF,1}=microbeID;
        gapfilledReactions{cntGF,2}='Enabling anaerobic growth';
        gapfilledReactions{cntGF,3}='Condition-specific gapfilling';
        gapfilledReactions(cntGF,4:length(condGF)+3)=condGF;
        cntGF=cntGF+1;
    end
    if ~isempty(targetGF)
        summary.targetGF=union(summary.targetGF,targetGF);

        gapfilledReactions{cntGF,1}=microbeID;
        gapfilledReactions{cntGF,2}='Enabling anaerobic growth';
        gapfilledReactions{cntGF,3}='Targeted gapfilling';
        gapfilledReactions(cntGF,4:length(targetGF)+3)=targetGF;
        cntGF=cntGF+1;
    end
    if ~isempty(relaxGF)
        summary.relaxGF=union(summary.relaxGF,relaxGF);

        gapfilledReactions{cntGF,1}=microbeID;
        gapfilledReactions{cntGF,2}='Enabling anaerobic growth';
        gapfilledReactions{cntGF,3}='Gapfilling based on relaxFBA';
        gapfilledReactions(cntGF,4:length(relaxGF)+3)=relaxGF;
        cntGF=cntGF+1;
    end
end

[AerobicGrowth, AnaerobicGrowth] = testGrowth(model, biomassReaction);
if AerobicGrowth(1,2) < tol
    % identify blocked reactions on complex medium
    model=useDiet(model,constraints);
    [model,condGF,targetGF,relaxGF] = runGapfillingFunctions(model,biomassReaction,biomassReaction,'max',database);
    % export the gapfilled reactions
    if ~isempty(condGF)
        summary.condGF=union(summary.condGF,condGF);

        gapfilledReactions{cntGF,1}=microbeID;
        gapfilledReactions{cntGF,2}='Growth on complex medium';
        gapfilledReactions{cntGF,3}='Condition-specific gapfilling';
        gapfilledReactions(cntGF,4:length(condGF)+3)=condGF;
        cntGF=cntGF+1;
    end
    if ~isempty(targetGF)
        summary.targetGF=union(summary.targetGF,targetGF);

        gapfilledReactions{cntGF,1}=microbeID;
        gapfilledReactions{cntGF,2}='Growth on complex medium';
        gapfilledReactions{cntGF,3}='Targeted gapfilling';
        gapfilledReactions(cntGF,4:length(targetGF)+3)=targetGF;
        cntGF=cntGF+1;
    end
    if ~isempty(relaxGF)
        summary.relaxGF=union(summary.relaxGF,relaxGF);

        gapfilledReactions{cntGF,1}=microbeID;
        gapfilledReactions{cntGF,2}='Growth on complex medium';
        gapfilledReactions{cntGF,3}='Gapfilling based on relaxFBA';
        gapfilledReactions(cntGF,4:length(relaxGF)+3)=relaxGF;
        cntGF=cntGF+1;
    end
end

% identify blocked biomass precursors on defined medium for the organism
[growsOnDefinedMedium,constrainedModel,~] = testGrowthOnDefinedMedia(model, microbeID, biomassReaction, inputDataFolder);
if growsOnDefinedMedium == 0
    % find reactions that are preventing the model from growing
    [model,condGF,targetGF,relaxGF] = runGapfillingFunctions(constrainedModel,biomassReaction,biomassReaction,'max',database,1);
    % export the gapfilled reactions
    if ~isempty(condGF)
        summary.condGF=union(summary.condGF,condGF);

        gapfilledReactions{cntGF,1}=microbeID;
        gapfilledReactions{cntGF,2}='Growth on defined medium';
        gapfilledReactions{cntGF,3}='Condition-specific gapfilling';
        gapfilledReactions(cntGF,4:length(condGF)+3)=condGF;
        cntGF=cntGF+1;
    end
    if ~isempty(targetGF)
        summary.targetGF=union(summary.targetGF,targetGF);

        gapfilledReactions{cntGF,1}=microbeID;
        gapfilledReactions{cntGF,2}='Growth on defined medium';
        gapfilledReactions{cntGF,3}='Targeted gapfilling';
        gapfilledReactions(cntGF,4:length(targetGF)+3)=targetGF;
        cntGF=cntGF+1;
    end
    if ~isempty(relaxGF)
        summary.relaxGF=union(summary.relaxGF,relaxGF);

        gapfilledReactions{cntGF,1}=microbeID;
        gapfilledReactions{cntGF,2}='Growth on defined medium';
        gapfilledReactions{cntGF,3}='Gapfilling based on relaxFBA';
        gapfilledReactions(cntGF,4:length(relaxGF)+3)=relaxGF;
        cntGF=cntGF+1;
    end
end

% if there are any false negative predictions
% find what reactions should be added to enable agreement with experimental
% data

fields=fieldnames(testResults);

for i=1:length(fields)
    % define if objective should be maximized or minimized
    if any(contains(fields{i},{'Carbon_sources','Metabolite_uptake','Drug_metabolism'}))
        osenseStr = 'min';
    elseif any(contains(fields{i},{'Fermentation_products','Secretion_products','Bile_acid_biosynthesis','PutrefactionPathways'}))
        osenseStr = 'max';
    end
    FNlist = testResults.(fields{i});
    if size(FNlist,1)>1
        FNs = FNlist(find(strcmp(FNlist(:,1),microbeID)),2:end);
        FNs = FNs(~cellfun(@isempty, FNs));
        if ~isempty(FNs)
            for j=1:length(FNs)
                metExch=['EX_' database.metabolites{find(strcmp(database.metabolites(:,2),FNs{j})),1} '(e)'];
                if contains(fields{i},'PutrefactionPathways')
                    % reaction ID itself provided
                    metExch = FNs{j};
                end
                % find reactions that could be gap-filled to enable flux
                try
                    [model,condGF,targetGF,relaxGF] = runGapfillingFunctions(model,metExch,biomassReaction,osenseStr,database);
                catch
                    condGF = {};
                    targetGF = {};
                    relaxGF = {};
                end
                % export the gapfilled reactions
                if ~isempty(condGF)
                    summary.condGF=union(summary.condGF,condGF);

                    gapfilledReactions{cntGF,1}=microbeID;
                    gapfilledReactions{cntGF,2}=FNs{j};
                    gapfilledReactions{cntGF,3}='Condition-specific gapfilling';
                    gapfilledReactions(cntGF,4:length(condGF)+3)=condGF;
                    cntGF=cntGF+1;
                end
                if ~isempty(targetGF)
                    summary.targetGF=union(summary.targetGF,targetGF);

                    gapfilledReactions{cntGF,1}=microbeID;
                    gapfilledReactions{cntGF,2}=FNs{j};
                    gapfilledReactions{cntGF,3}='Targeted gapfilling';
                    gapfilledReactions(cntGF,4:length(targetGF)+3)=targetGF;
                    cntGF=cntGF+1;
                end
                if ~isempty(relaxGF)
                    summary.relaxGF=union(summary.relaxGF,relaxGF);

                    gapfilledReactions{cntGF,1}=microbeID;
                    gapfilledReactions{cntGF,2}=FNs{j};
                    gapfilledReactions{cntGF,3}='Gapfilling based on relaxFBA';
                    gapfilledReactions(cntGF,4:length(relaxGF)+3)=relaxGF;
                    cntGF=cntGF+1;
                end
            end
        end
    end
end

% rebuild periplasmatic space if applies
infoFile = readInputTableForPipeline(infoFilePath);
[model] = createPeriplasmaticSpace(model,microbeID,infoFile);

% remove futile cycles if any exist
[model, deletedRxns, addedRxns] = removeFutileCycles(model, biomassReaction, database);
replacedReactions{1,1}=microbeID;
replacedReactions{1,2}='Futile cycle correction';
replacedReactions{1,3}='To replace';
replacedReactions(1,4:length(deletedRxns)+3)=deletedRxns;

% if any futile cycles remain
[atpFluxAerobic, atpFluxAnaerobic] = testATP(model);
if atpFluxAerobic > 200 || atpFluxAnaerobic > 150
    % let us try if running removeFutileCycles again will work
    [model, deletedRxns, addedRxns] = removeFutileCycles(model, biomassReaction, database);
    replacedReactions{1,1}=microbeID;
    replacedReactions{1,2}='Futile cycle correction';
    replacedReactions{1,3}='To replace';
    replacedReactions(1,4:length(deletedRxns)+3)=deletedRxns;
end

% add the gap-filling to model.comments field
for i=1:length(model.rxns)
    if strcmp(model.grRules{i},'demeterGapfill')
        model.grRules{i}=strrep(model.grRules{i},'demeterGapfill','');
        if ~isempty(find(strcmp(summary.condGF,model.rxns{i})))
            model.comments{i}='Added by DEMETER to enable flux with VMH-consistent constraints.';
        elseif ~isempty(find(strcmp(summary.targetGF,model.rxns{i})))
            model.comments{i}='Added by DEMETER during targeted gapfilling to enable production of required metabolites.';
        elseif ~isempty(find(strcmp(summary.relaxGF,model.rxns{i})))
            model.comments{i}='Added by DEMETER based on relaxFBA. Low confidence level.';
        end
    end
end

% rebuild and export the model
revisedModel = rebuildModel(model,database);
revisedModel = changeObjective(revisedModel,biomassReaction);
% limit flux through sink reactions
revisedModel = changeRxnBounds(revisedModel,revisedModel.rxns(find(strncmp(revisedModel.rxns,'sink_',5)),1),-1,'l');

end

author: opencobra

repo: cobratoolbox


text:
function plotlysetup_offline(plotly_bundle_url, varargin)

    % CALL: plotlysetup_offline(plotly_bundle_url);
    % WHERE: plotly_bundle_url is the plotly bundle url, e.g. http://cdn.plot.ly/plotly-latest.min.js
    % If no argument is provided, the default http://cdn.plot.ly/plotly-latest.min.js is used.
    % [1] adds plotly api to matlabroot/toolboxes. If successful do [2]
    % [2] adds plotly api to searchpath via startup.m of matlabroot and/or userpath

    %DEFAULT OUTPUT
    exception.message = '';
    exception.identifier = '';

    try %check number of inputs
        if(nargin == 0)
            plotly_bundle_url = 'http://cdn.plot.ly/plotly-latest.min.js';
        elseif (nargin>1)
            error('plotly:wrongInput',....
                ['\n\nWhoops! Wrong number of inputs. Please run >> help plotlysetup_offline \n',...
                'for more information regarding the setup your Plotly API MATLAB \n',...
                'Library. Please post a topic on https://community.plotly.com/c/api/matlab/ for more information.']);
        end
    catch exception %plotlysetup input problem catch...
        fprintf(['\n\n' exception.identifier exception.message '\n\n']);
        return
    end

    try
        %check to see if plotly is in the searchpath
        plotlysetupPath = which('plotlysetup');
        plotlyFolderPath = fullfile(fileparts(plotlysetupPath),'plotly');
        %if it was not found
        if (strcmp(genpath(plotlyFolderPath),''))
            error('plotly:notFound',...
                ['\n\nShoot! It looks like MATLAB is having trouble finding the current version '  ...
                '\nof Plotly. Please make sure that the Plotly API folder is in the same '  ...
                '\ndirectory as plotlysetup.m. Questions? Ask https://community.plotly.com/c/api/matlab/\n\n']);
        end
        %add Plotly API MATLAB Library to search path
        addpath(genpath(plotlyFolderPath));
    catch exception %plotly file not found problem catch
        fprintf(['\n\n' exception.identifier exception.message '\n']);
        return
    end

    if(~is_octave)

        try
            %embed the api to the matlabroot/toolbox dir.
            fprintf('\nAdding Plotly to MATLAB toolbox directory ...  ');

            %plotly folder in the matlab/toolbox dir.
            plotlyToolboxPath = fullfile(matlabroot,'toolbox','plotly');

            if(exist(plotlyToolboxPath,'dir')) %check for overwrite...
                fprintf(['\n\n[UPDATE]: \n\nHey! We see that a copy of Plotly has previously been added to\n' ...
                    'your Matlab toolboxes. Would you like us to overwrite it with:\n' plotlyFolderPath ' ? \n'...
                    'Careful! You may lose data saved to this Plotly directory.\n\n']);

                overwrite = input('Overwrite (y/n) ? : ','s');

                if(strcmpi(overwrite,'y'));
                    fprintf('\n[OVERWRITE]:\n\nOverwriting Plotly! ... Done \n');
                else
                    fprintf('\n[NO OVERWRITE]:\n\nDid not overwrite Plotly! ... Done \n');
                end
            else %toolbox Plotly not yet created

                %worked (without interuption)...just a formatting thing!
                fprintf('Done\n');

                %make the plotlyToolboxPath dir.
                status = mkdir(plotlyToolboxPath);

                %set status to overwrite
                overwrite = 'y';

                %check that the folder was created
                if (status == 0)
                    error('plotly:savePlotly', permissionMessage('save the Plotly folder'));
                end
            end

            if(strcmpi(overwrite,'y'))

                %move a copy of the Plotly api to matlab root directory
                [status, msg, messid] = copyfile(plotlyFolderPath,plotlyToolboxPath);
                %check that the plotly api was copied to the matlab root toolbox directory
                if (status == 0)
                    if(~strcmp(messid, 'MATLAB:COPYFILE:SourceAndDestinationSame'))
                        error('plotly:copyPlotly',permissionMessage('copy the Plotly folder'));
                    end
                end

            end

            %add it to the searchpath (startup.m will handle this next time!)
            addpath(genpath(plotlyToolboxPath),'-end');

            %save plotly api searchpath to startup.m files (only do this if we actually were able to store the api in mtlroot/toolbox!)
            fprintf('Saving Plotly to MATLAB search path via startup.m ... ');

            %check for a startup.m file in matlab rootpath (we want to add one here)
            startupFile = [];
            startupFileRootPath = fullfile(matlabroot,'toolbox','local','startup.m');
            if(~exist(startupFileRootPath,'file'))
                startFileID = fopen(startupFileRootPath, 'w');
                %startup.m does not exist and startupFilePath is non-writable
                if(startFileID == -1)
                    error('plotly:rootStartupCreation',permissionMessage('write the startup.m script'));
                end
                startupFile = {startupFileRootPath}; %needed because MATLAB only looks for startup.m when first opened.
            end

            %check for all startup.m file in searchpath
            startupFile = [startupFile; cell(which('startup.m','-all'))];
            %write the addpath - plotly api to the startup.m files
            [warnings] = addplotlystartup(startupFile);

            %worked!
            fprintf(' Done\n');

            %print any addplotlydstatup warnings;
            w = cellfun(@isempty,warnings);
            if(find(~w))
                %output warnings
                exception.warnings = warnings;
                fprintf(warnings{find(~w)});
            end

        catch exception %copying to toolbox/writing to startup.m permission problem catch...
            fprintf(['\n\n' exception.identifier exception.message '\n\n']);
        end

    else %if octave
        fprintf('\n\nOctave users: Automatic Plotly API embedding coming soon!\n\n');
    end %end check for matlab...

    %get offline bundle
    fprintf('\nNow downloading the plotly offline bundle ...');
    getplotlyoffline(plotly_bundle_url);

    %greet the people!
    fprintf('\nWelcome to Plotly! If you are new to Plotly please enter: >> plotlyhelp to get started!\n\n')

    end

    %helper message function
    function message = permissionMessage(spec)
    message = ['\n\nShoot! We tried to ' spec ' to the MATLAB toolbox \n',...
        'directory, but were denied write permission. You''ll have to add\n',...
        'the Plotly folder to your MATLAB path manually by running: \n\n',...
        '>> plotly_path = fullfile(pwd, ''plotly'')\n',...
        '>> addpath(genpath(plotly_path))\n\n',...
        'Questions? Ask https://community.plotly.com/c/api/matlab/\n\n'];
    end

author: plotly


text:
function [molStr,molTable] = getMolFormulaVMH(metabolite,metaboliteTable)
% This function returns the molecular formula of a given vmh metabolite.
%
% USAGE:
%       1) [molStr,molTable] = getMolFormulaVMH(metabolite)
%       2) [molStr,molTable] = getMolFormulaVMH(metabolite,metaboliteTable)
%
% INPUT:
%   metabolite:         a metabolite in the vmh database
%
% OPTIONAL INPUT:
%   metaboliteTable:    a metabolite table containing all vmh metabolites.
%                       This optional input is recommended if the function
%                       is run repetitively in a loop to eliminate
%                       readtable time.
%
% OUTPUT:
%   molStr:             a string containing the molecular formula of the
%                       metabolite
%   molTable:           a cellular array containing a breakdown of the
%                       frequency of each element in the molecule
%
% AUTHORS:
%       Bronson R. Weston, 2022

if ~exist('metaboliteTable','Var')
    metaboliteTable=readtable('MetaboliteDatabase.txt');
end
%try metabolite abreviation
ind=find(strcmp(metaboliteTable.Var1,metabolite));
if isempty(ind)
    %try metabolite full name
    ind=find(strcmp(metaboliteTable.Var2,metabolite));
    if isempty(ind)
        error('Invalid input metabolite!')
    elseif length(ind)>1
        ind=ind(1);
    end
end

%Get metabolite formula
molStr=metaboliteTable.Var4{ind};
uC=isstrprop(molStr,'upper');
tracker={};
ind=0;

%Break down metabolite by element and create molTable
for i=1:length(uC)
    if uC(i)==1
        ind=ind+1;
        tracker{ind}=molStr(i);

    else
        tracker{ind}=[tracker{ind},molStr(i)];
    end
end

molTable={};
for i=1:length(tracker)
    unit=tracker{i};
    num=isstrprop(unit,'digit');
    if ~any(num==1)
        molTable(i,:)={tracker{i} , 1};
    else
        molTable{i,1}=unit(num~=1);
        molTable{i,2}=str2num(unit(num==1));
    end
end
end


text:
%% Adaptative Quadrature
% This is an adaptative quadrature algorithm for Simpson and Trapezoidal
% rule.

clear
clc
close all

%% Quadrature grid and initial values
% Integration range and values
%
a=0; % Lower limit
b=1; % Upper limit
n=10; % Initial number of grid divisions
h=(b-a)/n; % Initial step size
x=zeros(1,n+1); % Quadrature points vector
epsilon=1/(b-a)*(1e-8); % Decided error

%% Stablish quadrature point x values
x(1)=a;
for i=1:n
    x(i+1)=x(i)+h;
end

%% Run quadrature function
for i=1:n
    S(i)=SimpsonAQ(x(i),x(i+1),epsilon); %Preseted function of AQ
end
total=sum(S)

repo: Aero-matlab