text:
% Simulates an MFSK system using the PDFs of the signals which are inputs to
% the envelope detector

%Proakis p.307
%Uses RS encoding and decoding

clear
clc

%***********************
%*** MFSK Parameters ***
%***********************
m = 4
M = 2^m         % Amount of Symbols in the MFSK system

E_s = 1        % Signal Power in dB

%**************************
%*** RS code Parameters ***
%**************************

%m = 8
n = 2^m - 1
k = 3                   %Information symbols [ 1 - n-1] where n = 2^m - 1
h = n-k
t = h/2

%*** Generate the Galois Field ***

field = gftuple([-1:2^m-2]', m, 2);

%*** Generator Polynomial ***

%Lin + Costello, p.171

%Get the generator polynomial
c = [1 0];
p(1) = c(1);

for i = 1:h-1
    p(1) = gfmul(p(1),1,field);
    p(2) = 0;
    c = gfconv(c,p,field);
end
g = c;

%**************************

%*****************************
%*** Simulation parameters ***
%*****************************

% number of iterations
num_runs = 25

% the length of the data to be transmitted over the channel
num_data = 100       %Number of codewords to be transmitted (Increase to 10000 to obtain smoother curves)

length_data_n = n * num_data;

%***************************

SNR_dB = 2
inc = 0.5;

T_SER = [];
ED_UC_SER = [];
ED_C_SER = [];

SNR_arr = [];

for nr = 1:num_runs

    nr

    SNR_dB = SNR_dB + inc;

    %Clear the counters
    ACC_Stats_ED = 0;
    ACC_Stats_ED_DEC = 0;

    % for each transmitted codeword, do the following

    for num_codewords = 1:num_data

        %Create k random symbols to be encoded
        INFO = randint(1,k,[-1 M-2]);

        %Encode the INFO to form a RS Codeword
        RS_CODE = RS_ENC(INFO,n,k,g,field);

        %convert data to M_FSK data
        send_MFSK = ConvertRS2MFSK(RS_CODE);

        %***********************************

        %***************************
        %*** Do the demodulation ***
        %***************************

        Metric_table = MFSK_DEMOD(send_MFSK, M,E_s, SNR_dB); %Create a table of metrics for MFSK with specific SNR

        RECEIVED_ED = [];

        for j = 1:n

            Metrics = transpose(Metric_table(1:M,j));

            %--- Do for Envelope Detection ---
            RECEIVED_ED = [RECEIVED_ED DECODE_MFSK(Metrics,M)];

        end

        %*** End demodulation ***

        %**********************
        %*** Do RS decoding ***
        %**********************

        %*** Envelope detection ***

        RECEIVED_ED_SYMB = ConvertMFSK2RS(RECEIVED_ED);
        Stats_ED =  Compare(RECEIVED_ED_SYMB, RS_CODE,n,field);       %Without coding

        Stats_ED_DEC = Decode_and_compare(RECEIVED_ED,RS_CODE,n,k,t,h,g,field);

        %**********************

        %************************
        %*** Accumulute Stats ***
        %************************
        ACC_Stats_ED = ACC_Stats_ED + Stats_ED;
        ACC_Stats_ED_DEC = ACC_Stats_ED_DEC + Stats_ED_DEC;

        %************************

    end

    %***********************
    %*** Calculate Stats ***
    %***********************

    %*** Theoretical result for envelope detection ***
    prob_symb_err = prob_symb_error_MFSK(M,SNR_dB);

    %*** Envelope Detector Statistics ***
    sym_prob_ED = ACC_Stats_ED / length_data_n;      %Envelope demodulation without decoding
    dec_sym_prob_ED = ACC_Stats_ED_DEC/length_data_n;    %Envelope demodulation with decoding

    %*******************

    T_SER = [T_SER prob_symb_err];
    ED_UC_SER = [ED_UC_SER sym_prob_ED];
    ED_C_SER = [ED_C_SER dec_sym_prob_ED];

    SNR_arr = [SNR_arr SNR_dB];

end

%Plot the data
plot_compare = figure;
semilogy(SNR_arr,T_SER,'k*-' , SNR_arr, ED_UC_SER,'kd-' , ...
    SNR_arr, ED_C_SER,'kx-');

xlabel('SNR','FontWeight','Bold','FontSize',12);
ylabel('Probability of a symbol error','FontWeight','Bold','FontSize',12);
tt = [int2str(2^m) '-FSK (' int2str(n) ',' int2str(k) ') RS']
title(tt);
set(plot_compare,'Position',[10,10,1000,600]);

legend('Theoretical','Envelope detection - uncoded', 'Envelope detection - coded')


text:
function X = Inv_SeparateAngles(R,deep,IsImageReal,MaxIts,ErrTol,X0)
% Inv_SeparateAngles: 'Invert' Separate Angles by the method of
%                           Least Squares
%  Usage:
%    X = Inv_SeparateAngles(R,deep,MaxIts,ErrTol);
%  Inputs:
%    R       Array of smoothly localized Fourier samples
%    deep    Depth of the angular splitting
%    MaxIts  Maximum number of CG iterations. Default 10.
%    ErrTol  Error tolerance. Default 1.e-9.
%  Outputs:
%    X   n by n matrix
%  Description
%    Performs the inverse angular separation. This is an
%    approximate inverse which uses a conjugate gradient solver on
%    the associated Gram system. The system is mildly preconditioned.
% See Also
%   SeparateAngles, Inv_AtA_CG,
%
% By Emmanuel Candes, 2003-2004

if nargin < 5,
	ErrTol = 1.e-9;
end
if nargin < 4,
	MaxIts = 10;
end
if nargin < 3,
  IsImageReal = 0;
end

n = size(R,3)*4;

W = SetScaleToZero(n);

P   = Adj_SeparateAngles(R,deep,IsImageReal).*W;
X   = Inv_AtA_CG(P,deep,MaxIts,ErrTol,X0);

author: thomaskuestner


text:
function toThing = copyfields(toThing, fromThing, varargin)
    % copy fields (or properties) from one structure (or object) to another. Types may be different
    %
    % toThing = copyfields(toThing, fromThing) copies fields/properties from fromThing into toThing,
    % returning toThing.  if THISONE is a handle object, the return can be ignored.
    %
    % Options:
    % copyfields(... , "safe") will copy all fields.  A warning will be displayed for any field that
    % doesn't exist.
    % copyfields(... , "safe", "quietfail") will not warn or error when a field doesn't exist.

    failQuietly =  any(varargin == "quietfail");

    if isstruct(fromThing)
        fn = fieldnames(fromThing);
    elseif isobject(fromThing)
        fn = properties(fromThing);
    elseif failQuietly
        return
    else
        error('cannot copy fields/properties from ["%s"], since it is neither an object nor a struct', inputname(2));
    end

    if any(varargin=="safe")
        for i = 1 : numel(fn)
            try
                toThing.(fn{i}) = fromThing.(fn{i});
            catch ME
                if failQuietly
                    do_nothing(); % do not copy
                else
                    warning(ME.identifier, 'Error copying property/field: %s',fn{i})
                end
            end
        end
    else
        for i = 1 : numel(fn)
            toThing.(fn{i}) = fromThing.(fn{i});
        end

    end
end

author: CelsoReyes


text:
function X_rec = recoverData(Z, U, K)
%RECOVERDATA Recovers an approximation of the original data when using the
%projected data
%   X_rec = RECOVERDATA(Z, U, K) recovers an approximation the
%   original data that has been reduced to K dimensions. It returns the
%   approximate reconstruction in X_rec.
%

% You need to return the following variables correctly.
X_rec = zeros(size(Z, 1), size(U, 1));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the approximation of the data by projecting back
%               onto the original space using the top K eigenvectors in U.
%
%               For the i-th example Z(i,:), the (approximate)
%               recovered data for dimension j is given as follows:
%                    v = Z(i, :)';
%                    recovered_j = v' * U(j, 1:K)';
%
%               Notice that U(j, 1:K) is a row vector.
%

X_rec = Z * U(:, 1:K)';

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

end

author: fewtime

repo: ML

save_path: github-repos/MATLAB/fewtime-ML


text:
%compute normalization of the spinelength by the mean spinelength size in the whole video
function [data,units]=compute_normlength(trx,n)
larvae= trx.exp2flies{n};
numlarvae = numel(larvae);

meanlength=zeros(1,numlarvae);
numelements=zeros(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    meanlength(i)=nanmean(trx(larva).spinelength(1:end));
    numelements(i)=numel(trx(larva).spinelength(1:end));
end
normvalue=sum(meanlength.*numelements./(sum(numelements)));
normlength=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    normlength{1,i}=trx(larva).spinelength(1:end)/normvalue;
end

units=parseunits([]);
data=normlength;

author: kristinbranson

repo: JAABA


text:
% Convert a stroke [x,y,t] such that it is uniformly sampled
% in time. This is done by linear interpolation in time
%
% Input
%  stk:  [n x 2] for x and y
%  time: [n x 1] for time coordinate
%  tint: [scalar] time interval (in milliseconds)
%
% Output
%  unif_stk: [k x 2] new stroke x and y
%  unif_t: [k x 1] uniform time interval
%
function [unif_stk,unif_t] = uniform_time_lerp(stk,time,tint)

    mint = min(time);
    maxt = max(time);

    unif_t = mint:tint:maxt;
    % make sure we don't leave out the last point
    if unif_t(end) ~= maxt
       unif_t = [unif_t maxt];
    end
    nt = length(unif_t);

    unif_stk = zeros(length(unif_t),2);
    for i=1:nt
        ti = unif_t(i);

        diff = time-ti;
        if any(diff==0)
            sel = stk(diff==0,:);
            unif_stk(i,:) = mean(sel,1);
            continue
        end

        % find the point before and after in time
        indx_gt = find(diff>0,1,'first');
        indx_lt = find(diff<0,1,'last');

        x_gt = stk(indx_gt,:);
        x_lt = stk(indx_lt,:);
        t_gt = time(indx_gt);
        t_lt = time(indx_lt);

        % Compute the linear interpolation
        frac = (ti - t_lt) ./ (t_gt - t_lt);
        assert(frac<=1);
        unif_stk(i,:) = (1-frac).*x_lt + frac.*x_gt;
    end

end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
%% housekeeping
close all
clearvars
clc

%% read the model
m=rise('kaiji1207_exp');

%% alternatively, solve the model directly

m=solve(m);

%% print the solution of the model

m.print_solution
% alternative call
% print_solution(kaiji)

%% print solution of a subset of variables

m.print_solution({'C','H','K'})

%% doing things the RBC way: waste of time
% RISE write the solution of the whole system as: X_t=T*X_{t-1}+R*e_t. The
% RBC guys write solutions as S_t=P*S_{t-1}+K*e_t and Y_t=F*S_t . Below, we
% show how to recover matrices P, K and F.
% collect the T and R matrices
[T,R]=load_solution(m,'iov');
T=T{1}; R=R{1};
% set some elements to 0 sharp
T(abs(T)<1e-10)=0;
R(abs(R)<1e-10)=0;
% get the state columns
state_cols=any(T);
% Kaiji wastes time double-checking check that find(any(T)) gives the
% columns of the states
% get the control columns as the columns that are not states
control_cols=~state_cols;
% now build the P matrix expressing the states as a function of the states
P=T(state_cols,state_cols);
% recover the solution of the controls as a function of states
FP=T(control_cols,state_cols);
F=FP/P; % F=FP*inv(P) P should always be invertible unless collinearity
% recover the shock impact on the states
K=R(state_cols,:);

allvars=m.endogenous.name;

state_variables=allvars(state_cols);

control_variables=allvars(control_cols);

%% compute impulse responses
% In general, if you want to see the options for a particular method,
% create an empty rise object. e.g. tmp=rise.empty(0); then you can go
% ahead and call the specific function you want to use on that empty object
% and it will list out the various options and the defaults of those
% options. e.g. irf(tmp)
simple_irfs=irf(m,'irf_periods',20);
%% construct a vector of models with different anticipation options
kaiji_unant=m.set('irf_anticipate',false);
myvector=[m,kaiji_unant];

%% change some option in the vector
% here we set:
% 1- the number of irf periods. If we don't the default of 40 will be used
% 2- we assume that we know as of today that shocks will hit 3 periods from
% now. And we choose to act on that information (anticipated) or disregard
% the information (unanticipated)
% 3- we can also choose the sign of the shocks. By default, the shocks will
% be positive.
myirfs=irf(myvector,'irf_periods',20,'solve_shock_horizon',3,'irf_shock_sign',1);

%% plot the irfs
shock_list=m.exogenous.name;
tex=get(m,'tex');
% var_list=allvars;
% just re-ordering the variables in a way that I like for the plotting
var_list={'A','B','D','PSI','C','K','H','R'};
for ishock=1:numel(shock_list)
    shock=shock_list{ishock};
    figure('name',['IRFs to a ',tex.(shock), 'shock']);
    for ivar=1:numel(var_list)
        endovar=var_list{ivar};
        subplot(3,3,ivar)
        plot(myirfs.(shock).(endovar),'linewidth',2)
        title(tex.(endovar))
        if ivar==1
            legend('anticipated','unanticipated')
        end
    end
end
%% estimation. Now the DSGE way

%% read the data and create the time series
[datta,names]=xlsread('data.xlsx');
names=names(1,2:end);
start_date='1960';

mydata=struct();
for ii=1:numel(names)
    mydata.(names{ii})=ts(start_date,datta(:,ii+1));
end

%% plot the data
figure('name','Observed data');
for ii=1:numel(names)
    subplot(3,1,ii)
    plot(mydata.(names{ii}))
    title(names{ii});
end

%% estimate the model

m=estimate(m,'data',mydata);

%% historical decomposition of shocks
histdec=historical_decomposition(m);

%% plot the decomposition
figure('name','historical decomposition of shocks and initial conditions')
for ivar=1:numel(var_list)
    vname=var_list{ivar};
    subplot(3,3,ivar)
    plot_decomp(histdec.(vname))
    title(tex.(vname))
    if ivar==1
        contrib_names=histdec.(vname).varnames;
        for jj=1:numel(contrib_names)
            if ~isfield(tex,contrib_names{jj})
                continue
            end
            contrib_names{jj}=tex.(contrib_names{jj});
        end
        hleg=legend(contrib_names,...
            'Location','BestOutside','orientation','horizontal');
        pp=get(hleg,'position');
        pp(1:2)=0;
        set(hleg,'position',pp)
    end
end

%% counterfactual: what if only one shock had been alive?
for ishock=1:numel(shock_list)
    %,{'EPS_PSI','EPS_B','EPS_D'}
    [counterf,actual]=counterfactual(m,[],1,shock_list{ishock});
    figure('name',['Counterfactual: ',tex.(shock_list{ishock}),' shock only'])
    for ivar=1:numel(var_list)
        vname=var_list{ivar};
        subplot(3,3,ivar)
        plot([actual.(vname),counterf.(vname)])
        title(tex.(vname))
        if ivar==1
            legend({'actual','counterfactual'})
        end
    end
end

%% counterfactual for a subset of shocks
close all
[counterf,actual]=counterfactual(m,[],1,{'EPS_PSI','EPS_B','EPS_D'});
figure('name','Counterfactual: EPS_PSI and EPS_B shock only');
for ivar=1:numel(var_list)
    vname=var_list{ivar};
    subplot(3,3,ivar)
    plot([actual.(vname),counterf.(vname)])
    title(tex.(vname))
    if ivar==1
        legend({'actual','counterfactual'})
    end
end

%% variance decomposition
vardec=variance_decomposition(m);

%% plot the decomposition
close all
figure('name','Variance decomposition of shocks')
for ivar=1:numel(var_list)
    vname=var_list{ivar};
    subplot(3,3,ivar)
    plot_decomp('0:50',vardec.conditional.(vname))
    title(tex.(vname))
    if ivar==1
        contrib_names=vardec.conditional.(vname).varnames;
        for jj=1:numel(contrib_names)
            if ~isfield(tex,contrib_names{jj})
                continue
            end
            contrib_names{jj}=tex.(contrib_names{jj});
        end
        hleg=legend(contrib_names,...
            'Location','BestOutside','orientation','horizontal');
        pp=get(hleg,'position');
        pp(1:2)=0;
        set(hleg,'position',pp)
    end
end

author: jmaih


text:
function Xr = pcaRandVec( U, mu, vars, k, n, hypershpere, show )
% Generate random vectors in PCA subspace.
%
% Used to generate random vectors from the subspace spanned by the first k
% principal components.  The points generated come from the gaussian
% distribution from within the subspace. Can optionally generate points on
% the subspace that are also on a hypershpere centered on the origin.  This
% may be useful if the original data points were all from a hypershpere --
% for example they were normalized via imNormalize.  Set the optional
% hypershpere flag to 1 to generate points only on the hypersphere.
%
% USAGE
%  Xr = pcaRandVec( U, mu, vars, k, n, [hypershpere], [show] )
%
% INPUTS
%  U           - returned by pca.m
%  mu          - returned by pca.m
%  vars        - returned by pca.m
%  k           - number of principal coordinates to use
%  n           - number of points to generate
%  hypershpere - [0] generate points on hypersphere (see above)
%  show        - [1] figure to use for display (no display if == 0)
%
% OUTPUTS
%  Xr          - resulting randomly generated vectors
%
% EXAMPLE
%
% See also PCA IMNORMALIZE
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<6 || isempty(hypershpere) ); hypershpere=0; end
if( nargin<7 || isempty(show) ); show=0; end

siz1 = size(mu);   nd = ndims(mu);
sizX = [siz1, n];  D = prod(siz1);

% generate random vectors inside of subspace.
C = diag( vars(1:k).^-1 );
Yr = C * randn(k,n);
Uk = U(:,1:k);
Xr = Uk * Yr;

if(hypershpere)
  % The final point must lie on hypershpere.  Need to scale each element xr
  % of Xr so that after adding it to mu the resulting vector has length D.
  % So need to find k>0 such that (k*xr + mu) has length D.  To find k
  % simply solve the quadratic equation induced by ||(k*xr + mu)||=D,
  % choosing the root such that k>0.  Note that the mean of resulting
  % vector will not necessarily be 0, but the variance will pretty much be
  % 1! Regardless, renomralize at the end.
  muv = mu(:);  muMag = dot(muv,muv);
  for i=1:n
    xr = Xr(:,i);
    rs = roots( [dot(xr,xr), 2 * dot(xr,muv), muMag-D] );
    Xr(:,i) = muv + max(rs)*xr;
  end
  Xr = reshape( Xr, sizX );
  Xr = fevalArrays( Xr, @imNormalize );
else
  % simply add the mean to reshaped Xr
  Xr = reshape( Xr, sizX );
  muRep = repmat(mu, [ones(1,nd), n ] );
  Xr = Xr + muRep;
end

% optionaly show resulting vectors
if(show && (nd==2 || nd==3)); figure(show); montage2(Xr); end

%%% Little test - see if eigenvectors induced by randomly generated vectors
%%% are the same as the original eigenvectors.  Assumes [U,mu,vars] exist.
%   Xr = pcaRandVec( U, mu, vars, 3, 100 );
%   [ Ur, mur, varsr ] = pca( Xr );
%   ind = 3;
%   Uim = reshape( U(:,ind), [ size(mu,1), size(mu,2) ]  );
%   Uimr = reshape( Ur(:,ind), [ size(mu,1), size(mu,2) ]  );
%   if( sum(abs(Uim-Uimr))>sum(abs(Uim+Uimr))) Uimr=Uimr*-1; end; %sign?
%   clf; subplot(3,1,1); im( Uim ); subplot(3,1,2); im( Uimr );
%   subplot(3,1,3); im( Uim - Uimr);
%%%

author: zouchuhang

repo: LayoutNet


text:
%PRCURSOR Show object ident.
%
%    PRCURSOR(H)
%
% Enable the datacursor in a scatterplot. This can be used to
% investigate the object identifier by clicking on the object.

% Copyright: D.M.J. Tax, D.M.J.Tax@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function out = prcursor(h,event_obj)

% This file has to tasks:
% 1. to setup the datacursor (by enabling it and setting the callback
% function)
% 2. to provide the callback when a user pressed on an object

if ~exist('datacursormode','file')
	error('MATLAB 7.0 or newer is required (datacursormode.m is not available).');
end
if nargin ~=2
	% we are doing the setup:
	if nargin<1
		h = gcf;
	end
	dh = datacursormode(h);
	set(dh,'enable','on','updatefcn',@prcursor);

else
	% we are doing the callback:
	nr = get(event_obj,'dataindex');
	ud = get(get(event_obj,'target'),'UserData');
	if ~isempty(ud) & isfield(ud,'ident')
		nr = ud.ident(nr);
	end
	out = sprintf('obj. %d',nr);
end

return

author: marianux

repo: ecg-kit

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


text:
function G=doog2(sig,r,th,N);
% G=doog2(sig,r,th,N);
% Make difference of offset gaussians kernel
% theta is in degrees
% (see Malik & Perona, J. Opt. Soc. Amer., 1990)
%
% Example:
% >> imagesc(doog2(1,12,0,64,1))
% >> colormap(gray)
%
% Serge Belongie

no_pts=N;  % no. of points in x,y grid

[x,y]=meshgrid(-(N/2)+1/2:(N/2)-1/2,-(N/2)+1/2:(N/2)-1/2);

phi=pi*th/180;
sigy=sig;
sigx=r*sig;
R=[cos(phi) -sin(phi); sin(phi) cos(phi)];
C=R*diag([sigx,sigy])*R';

X=[x(:) y(:)];

Gb=gaussian(X,[0 0]',C);
Gb=reshape(Gb,N,N);

m=R*[0 sig]';
Ga=gaussian(X,m,C);
Ga=reshape(Ga,N,N);
Gc=rot90(Ga,2);

a=-1;
b=2;
c=-1;

G = a*Ga + b*Gb + c*Gc;

author: jwyang


text:
function check_FourierGrad

S = mtexdata('dubna');
pf = S.pf;
odf = calcODF(pf,'zeroRange');
fodf = FourierODF(odf);

ori = orientation.rand(1000,pf.CS);
g1 = fodf.grad(ori(:),'check','delta',0.1*degree);
g2 = fodf.grad(ori(:));

if max(norm(g1-g2)./norm(g1)) < 1e-2
  disp(' Fourier gradient test passed');
else
  disp(' Fourier gradient test failed');
end

return

%% test 1 - check f_theta
cs = crystalSymmetry('1');
odf = fibreODF(Miller(0,0,1,cs),vector3d.Z,'halfwidth',40*degree);
fodf = FourierODF(odf);

omega = linspace(-179,179) * degree;
ori = rotation.byAxisAngle(vector3d(0,1,0),omega);

g1 = fodf.grad(ori(:),'check','delta',0.1*degree);
g2 = fodf.grad(ori(:));

if max(norm(g1-g2)) < 1e-3, disp(' Test 1 passed'); end

%% test 2

ori = rotation.byAxisAngle(vector3d(1,2,0)s,omega);

g1 = fodf.grad(ori(:),'check','delta',0.1*degree);
g2 = fodf.grad(ori(:));

if max(norm(g1-g2)) < 1e-3, disp(' Test 2 passed'); end

%% test 3 - check last line in grad composition

ori = rotation.byAxisAngle(vector3d(1,2,3),omega);

g1 = fodf.grad(ori(:),'check','delta',0.1*degree);
g2 = fodf.grad(ori(:));

if max(norm(g1-g2)) < 1e-2
  disp(' Test 3 passed');
else
  disp(' Test 3 failed');
end

%% test 4 - final f_theta test

ori = orientation.rand(100,cs);
g1 = fodf.grad(ori(:),'check','delta',0.1*degree);
g2 = fodf.grad(ori(:));

if max(norm(g1-g2)) < 1e-2
  disp(' Test 4 passed');
else
  disp(' Test 4 failed');
end

%% test 5 - f_phi1 test

odf = fibreODF(Miller(0,1,0,cs),vector3d.Y,'halfwidth',40*degree);
fodf = FourierODF(odf);

ori = orientation.byEuler(omega,89.9999*degree,0*degree,cs,'ABG');

%odf.grad(rot)
g1 = fodf.grad(ori(:),'check','delta',0.1*degree);
g2 = fodf.grad(ori(:));

if max(norm(g1-g2)) < 1e-2
  disp(' f_phi1 test passed');
else
  disp(' f_phi1 test failed');
end

%% test 6 - f_phi1 test

odf = fibreODF(Miller(0,1,0,cs),vector3d.Y,'halfwidth',40*degree);
fodf = FourierODF(odf);

ori = orientation.byEuler(0,89.9999*degree,omega,cs,'ABG');

%odf.grad(rot)
g1 = fodf.grad(ori(:),'check','delta',0.1*degree);
g2 = fodf.grad(ori(:));

if max(norm(g1-g2)) < 1e-2
  disp(' f_phi2 test passed');
else
  disp(' f_phi2 test failed');
end

%% test 5 - f_phi1 + f_phi_2 test

odf = fibreODF(Miller(0,1,0,cs),vector3d.Y,'halfwidth',40*degree);
fodf = FourierODF(odf);

ori = orientation.byAxisAngle(vector3d(0,1,100),omega);

%odf.grad(rot)
g1 = fodf.grad(ori(:),'check','delta',0.1*degree);
g2 = fodf.grad(ori(:));

if max(norm(g1-g2)) < 1e-2
  disp(' Test 6 passed');
else
  disp(' Test 6 failed');
end

%% test 6 - f_phi1 - f_phi_2 test

odf = fibreODF(Miller(0,1,0,cs),vector3d.Y,'halfwidth',40*degree);
fodf = FourierODF(odf);

ori = orientation.byAxisAngle(vector3d(0,1,100),omega);

%odf.grad(rot)
g1 = fodf.grad(ori(:),'check','delta',0.1*degree);
g2 = fodf.grad(ori(:));

if max(norm(g1-g2)) < 1e-2
  disp(' Test 6 passed');
else
  disp(' Test 6 failed');
end

%%
function testing

cs = crystalSymmetry('1');
%ref = orientation.id(cs);
ref = orientation.rand(cs);
odf = unimodalODF(ref,'halfwidth',40*degree);

fodf = FourierODF(odf)

omega = linspace(0,40)*degree;
%omega = 5*degree;
rot = ref*rotation(vector3d(1,1,1),omega)

g2 = fodf.grad(rot,'check','delta',0.05*degree)
g3 = odf.grad(rot)
g1 = fodf.grad(rot)

fodf.eval(rot)
odf.eval(rot)

%plot(norm(odf.grad(rot))

end

%%
function test2

cs = crystalSymmetry('1');
odf = fibreODF(Miller(0,1,0,cs),vector3d.Y,'halfwidth',40*degree);
fodf = FourierODF(odf)

omega = 50*degree;
%omega = linspace(0,179)*degree;
%ori = orientation('axis',vector3d(0,1,100),'angle',omega,cs);
%ori = orientation.byEuler(0*degree,0.01*degree,omega,cs,'ABG');

% f_phi1 test with cos(Phi) == 0 --> ok
ori = orientation.byEuler(omega,89.9999*degree,30*degree,cs,'ABG')

%%
% f_phi2 test with cos(Phi) == 1

omega = 80*degree
%omega = linspace(0,89)*degree;

odf = fibreODF(Miller(0,1,0,cs),vector3d.Z,'halfwidth',40*degree);
fodf = FourierODF(odf)

%ori = orientation.byEuler(0*degree,omega,[ 90]*degree,cs,'ABG');

%ori = orientation.rand(100,cs);

ori = orientation.byEuler(0*degree,90*degree,45*degree,cs,'ABG');
%ori = orientation.byEuler(10*degree,80*degree,[ 90]*degree,cs,'ABG');

g2 = fodf.grad(ori(:))
g1 = odf.grad(ori(:))

d_Phi(odf,ori)
d_phi1(odf,ori)
d_phi2(odf,ori)

%%

fhat1 = zeros(3);% fhat1([4 6]) = [1 -1]*1i;
fhat1 = zeros(3); %fhat1(9) = 1;
fhat2 = zeros(5); %fhat2([6,8,10]) = 1;
fhat2 = zeros(5); fhat2(4) = 1;
fhat3 = zeros(7); %fhat3([22 24 26 28]) = [1 1 -1 -1]*1i;
odf2 = FourierODF([1;fhat1(:);fhat2(:);fhat3(:)],cs);

[g1,d1,d2,d3] = odf2.components{1}.grad(ori3);
g1
g2 = odf2.grad(ori3,'check')

[d_phi1(odf2,ori3),d_Phi(odf2,ori3),d_phi2(odf2,ori3);...
  d1,d2,d3]

%%

%plot([norm(g1-g2)])
plot([g1.z , g2.z]) % x anfangs ok dann zu schnell negativ
plot([g1.y,g2.y])   % y - sollte konstant 0 sein
plot([g2.z,g1.z])  % z - faktor 4 zu groß
fodf.eval(ori)
odf.eval(ori)

%plot(norm(odf.grad(rot))

end

function test1

cs = crystalSymmetry('1');
odf = fibreODF(Miller(0,0,1,cs),vector3d.Z,'halfwidth',40*degree);
fodf = FourierODF(odf)

omega = 55*degree
%omega = linspace(-179,179) * degree;
ori = rotation.byAxisAngle(vector3d(-1,1,3),omega)
%ori = orientation.rand(cs);

%odf.grad(rot)
g1 = fodf.grad(ori(:),'check','delta',0.1*degree)
g2 = fodf.grad(ori(:))

plot([g1.y,g2.y])

end

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function vl_demo_gmm_3d
% VL_DEMO_GMM_3D   Show how to initialize GMM learning with a custom set of modes

%% Initialize points clustered in 3D
numClusters = 5 ;
dimension = 3 ;
cc = hsv(numClusters) ;
randn('state',0) ;
X = [] ;
for i=1:numClusters
  delta = 5*[i;0;0] ;
  %num = ceil(1000 / i) ;
  num = 500 ;
  s = diag([1,i, 1]) ;
  X = [X, bsxfun(@plus, s*randn(dimension, num), delta)] ;
end

%% Cluster using KMeans (ANN algorithm)
elapsedKMEANS = tic ;
vl_twister('state',0) ;
[initMeans, assignments] = vl_kmeans(X, numClusters, ...
                                     'Verbose', ...
                                     'Algorithm', 'ann', ...
                                     'MaxNumIterations', 10, ...
                                     'Distance', 'l2', ...
                                     'NumTrees', 3, ...
                                     'MaxNumComparisons', 10);
elapsedKMEANS = toc(elapsedKMEANS) ;
fprintf('\n ---------- \n KMeans time - %f \n ---------- \n\n',elapsedKMEANS) ;

%% Create initial GMM modes from KMeans clusters
initSigmas = zeros(dimension,numClusters);
initWeights = zeros(1,numClusters);

figure(1) ; clf ;
subplot(1,2,1) ; hold on ;

for i=1:numClusters
  Xk = X(:,assignments==i) ;
  initWeights(i) = size(Xk,2) / numClusters ;

  plot3(Xk(1,:),Xk(2,:),Xk(3,:),'.','color',cc(i,:)) ;
  if size(Xk,1) == 0 || size(Xk,2) == 0
    initSigmas(:,i) = diag(cov(X'));
  else
    initSigmas(:,i) = diag(cov(Xk'));
  end
end

axis equal ;
xlabel('x') ; ylabel('y') ; zlabel('z') ;
title('GMM: KMeans intialization') ;

%% Fit a GMM model to the maximum likelihood
elapsedGMM = tic ;
[means,sigmas,weights,ll,posteriors] = vl_gmm(X, numClusters, ...
                                              'Initialization','custom', ...
                                              'InitMeans',initMeans, ...
                                              'InitCovariances',initSigmas, ...
                                              'InitPriors',initWeights, ...
                                              'Verbose', ...
                                              'MaxNumIterations', 100) ;
elapsedGMM = toc(elapsedGMM);
fprintf('\n ---------- \n GMM time - %f \n ---------- \n',elapsedGMM) ;

subplot(1,2,2) ; hold on ;
[~,idx] = max(posteriors,[],1);
for i=1:numClusters
  plot3(X(1,idx == i),X(2,idx == i),X(3,idx == i),'.','color',cc(i,:)) ;
end

axis equal ;
xlabel('x') ; ylabel('y') ; zlabel('z') ;
title('GMM: Estimated Gaussian mixture') ;

author: jianxiongxiao

repo: ProfXkit


text:
function str = finish(str,unknowns)
%TEMPLATE/FINISH Apply given strategy to unknown fields in a string
%  STR = FINISH(STR,UNKNOWNS) applies on string STR the strategy defined
%  in UNKNOWNS to unknowns fields '{UNKNOWNS_FIELDS}'.
%  UNKNOWNS may be:
%    * 'keep' to do nothing
%    * 'remove' to remove all undefined fields
%    * 'comment' to replace undefined fields by a warning HTML comment.
%  This function uses Matlab REGEXPREP function coming with R13. If you
%  hold an older version, please comment lines 38 and 42: then you can
%  only apply the 'keep' strategy.

%  Copyright (C) 2003 Guillaume Flandin <Guillaume@artefact.tk>
%  $Revision: 1.0 $Date: 2003/05/05 22:19:51 $

error(nargchk(2,2,nargin));

switch lower(unknowns)
	case 'keep'
		%- do nothing
	case 'remove'
		%%%%%%%%%%%%%%%%%%%%%%%% WIH REGEXP ONLY %%%%%%%%%%%%%%%%%%%%
		% str = regexprep(str,'{[^ \t\r\n}]+}','');
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	case 'comment'
		%%%%%%%%%%%%%%%%%%%%%%%% WIH REGEXP ONLY %%%%%%%%%%%%%%%%%%%%
		% str = regexprep(str,'{[^ \t\r\n}]+}','<!-- Template variable undefined -->');
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	otherwise
		error('[Template] Unknown action.');
end


text:
function spec_plot2(varargin);
%
% Type: spec_plot2(S_m,freq,time,typ,f_min,f_max,c_min,c_max);
% Type: spec_plot2(S_m,freq,time,typ,f_min,f_max,c_min);
% Type: spec_plot2(S_m,freq,time,typ,f_min,f_max);
% Type: spec_plot2(S_m,freq,time,typ,f_min);
% Type: spec_plot2(S_m,freq,time,typ);
% Type: spec_plot2(S_m,freq,time);
%
% Inputs:
%
% S_m    := Time frequency distribution matrix n x m matrix
% freq   := Frequency vector n x 1 vector
% time   := Time vector m x 1 vector
% typ    := 'lin' or 'log' string for color scale of image plot
% f_min  := Minimum frequency to show, scalar
% f_max  := Maximum frequency to show, scalar
% c_min  := value for low end of color spectrum (see caxis), scalar
% c_max  := value for high end of color spectrum (see caxis), scalar
%
% This version uses pcolor and shading interp instead of imagesc

% Scot McNeill, University of Houston, Fall 2007.
%
msg=nargchk(3,8,nargin);
if ~isempty(msg)
   error(msg)
end
%
S_m=varargin{1};
freq=varargin{2};
time=varargin{3};
%
if (length(varargin) >= 4 & ~isempty(varargin{4}))
   typ=varargin{4};
else
   typ='log';
end
%
if (length(varargin) >= 5 & ~isempty(varargin{5}))
   f_min=varargin{5};
else
   f_min=min(freq);
end
%
if (length(varargin) >=6 & ~isempty(varargin{6}))
   f_max=varargin{6};
else
   f_max=max(freq);
end
%
if (length(varargin) >= 7 & ~isempty(varargin{7}))
   c_min=varargin{7};
else
   c_min=[];
end
%
if (length(varargin) ==8 & ~isempty(varargin{8}))
   c_max=varargin{8};
else
   c_max=[];
end
%
[nr,nc]=size(S_m);
nt=length(time);
nf=length(freq);
if (nr ~= nf & nc ~= nt)
   error('S_m must have dimensions length(freq) x length(time)')
end
%
[i1]=find(freq>=f_min & freq<=f_max);
freq=freq(i1);
S_m=S_m(i1,:);
%
newplot;
if strcmp(typ,'log');
   %imagesc(time,freq,log10(abs(S_m)+eps));
   pcolor(time,freq,log10(abs(S_m)+eps));shading interp;
   axis xy;
   colormap(jet);
   xlabel('Time (seconds)');
   ylabel('Frequency (Hz)');
   title('Spectrogram');
   axis([min(time),max(time),f_min,f_max]);
   [c1,c2]=caxis;
   if isempty(c_min);c_min=c1;end;
   if isempty(c_max);c_max=c2;end;
   caxis([c_min,c_max]);
   h = colorbar;
   set(get(h,'YLabel'),'String','Spectrogram, log_1_0([Units]^2/Hz)');
elseif strcmp(typ,'lin');
   %imagesc(time,freq,abs(S_m));
   pcolor(time,freq,abs(S_m));shading interp;
   axis xy;
   colormap(jet);
   xlabel('Time (seconds)');
   ylabel('Frequency (Hz)');
   title('Spectrogram');
   axis([min(time),max(time),f_min,f_max]);
   [c1,c2]=caxis;
   if isempty(c_min);c_min=c1;end;
   if isempty(c_max);c_max=c2;end;
   caxis([c_min,c_max]);
   h = colorbar;
   set(get(h,'YLabel'),'String','Spectrogram, [Units]^2/Hz');
else
   error(['typ must be ''lin'' or ''log''.']);
end
set(gca,'tickdir','out');
drawnow

author: Sable

repo: mcbench-benchmarks


text:
function context = document_context(name, varargin)
% document_context Create report context structure
%
% Creates a new report context structure that can be used to generate
% report documents.
%
% Input:
% - name (string): Name of the report. Is used to name a directory where the report is stored in.
% - varargin[Title] (string): A human-friendly name of the report. May be more verbose than report name.
% - varargin[FigureEPS] (boolean): Should the report also generate formaps in Encapsulated PostScript (for LaTeX).
% - varargin[FigureRaw] (boolean): Should internal raw figures and data also be saved for future processing.
% - varargin[Standalone] (boolean): If true (default), the JavaScript and CSS resources are copied to the report.
% - varargin[Cache] (boolean): Should cache be used.
%
% Output:
% - context (structure): A report context structure.
%

title = name;
figure_eps = false;
figure_raw = false;
standalone = true;
cache = false;

for i = 1:2:length(varargin)
    switch lower(varargin{i})
        case 'title'
            title = varargin{i+1};
        case 'figureeps'
            figure_eps = varargin{i+1};
        case 'figureraw'
            figure_raw = varargin{i+1};
        case 'standalone'
            standalone = varargin{i+1};
        case 'cache'
            cache = varargin{i+1};
        otherwise
            error(['Unknown switch ', varargin{i}, '!']) ;
    end
end

context.root = fullfile(get_global_variable('directory'), 'reports', name);
context.images = fullfile(context.root, 'images');
context.raw = fullfile(context.root, 'raw');
context.data = fullfile(context.root, 'data');
context.exporteps = figure_eps;
context.exportraw = figure_raw;
context.standalone = standalone;
context.prefix = '';
context.imagesurl = 'images';
context.rawurl = 'raw';
context.dataurl = 'data';
context.title = title;
context.cache = cache;

mkpath(context.root);
mkpath(context.images);
mkpath(context.data);
mkpath(context.raw);

context.cachedir = fullfile(context.root, 'cache');
mkpath(context.cachedir);

author: votchallenge

repo: toolkit-legacy


text:
function [f,J] = spm_fx_lfp(x,u,P,M)
% state equations for a neural mass model of erps
% FORMAT [f,J] = spm_fx_lfp(x,u,P,M)
% x      - state vector
%   x(:,1)  - voltage (spiny stellate cells)
%   x(:,2)  - voltage (pyramidal cells)         +ve
%   x(:,3)  - voltage (pyramidal cells)         -ve
%   x(:,4)  - current (spiny stellate cells)    +ve
%   x(:,5)  - current (pyramidal cells)         +ve
%   x(:,6)  - current (pyramidal cells)         -ve
%   x(:,7)  - voltage (inhibitory interneurons) +ve
%   x(:,8)  - current (inhibitory interneurons) +ve
%   x(:,9)  - voltage (pyramidal cells)
%   x(:,10) - voltage (inhibitory interneurons) -ve
%   x(:,11) - current (inhibitory interneurons) -ve
%   x(:,12) - voltage (inhibitory interneurons)
%
%   x(:,13) - slow potassium conductance
%
% f    = dx(t)/dt  = f(x(t))
% J    = df/dx
%
% Fixed parameter scaling [Defaults]
%
%  E = [32 16 4];             % extrinsic rates (forward, backward, lateral)
%  G = [1 1 1/2 1/2 1/8]*128; % intrinsic rates (g1, g2, g3, g4, g5)
%  D = [2 16];                % propagation delays (intrinsic, extrinsic)
%  H = [4 32];                % receptor densities (excitatory, inhibitory)
%  T = [4 16];                % synaptic constants (excitatory, inhibitory)
%  R = [2 1];                 % parameters of static nonlinearity
%
%__________________________________________________________________________
%
% This is a simplified version of spm_fx_erp
%
% David O, Friston KJ (2003) A neural mass model for MEG/EEG: coupling and
% neuronal dynamics. NeuroImage 20: 1743-1755
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_fx_lfp.m 5369 2013-03-28 20:09:27Z karl $

% check if intrinsic connections are free parameters
%--------------------------------------------------------------------------
try, P.H; catch, P.H = 0; end

% get dimensions and configure state variables
%--------------------------------------------------------------------------
x    = spm_unvec(x,M.x);       % neuronal states
n    = size(x,1);              % number of sources
s    = size(x,2);              % number of states

% [default] fixed parameters
%--------------------------------------------------------------------------
E    = [32 16 4];              % extrinsic rates (forward, backward, lateral)
G    = [1 1 1/2 1/2 1/32]*128; % intrinsic rates (g1, g2 g3, g4)
D    = [2 4];                  % propagation delays (intrinsic, extrinsic)
H    = [8 32];                 % receptor densities (excitatory, inhibitory)
T    = [4 16];                 % synaptic constants (excitatory, inhibitory)
R    = [1 2];                  % parameters of static nonlinearity

% [specified] fixed parameters
%--------------------------------------------------------------------------
if isfield(M,'pF')
    try, E  = M.pF.E; end
    try, G  = M.pF.H; end
    try, D  = M.pF.D; end
    try, H  = M.pF.G; end
    try, T  = M.pF.T; end
    try, R  = M.pF.R; end
end

% exponential transform to ensure positivity constraints
%--------------------------------------------------------------------------
A{1} = exp(P.A{1})*E(1);
A{2} = exp(P.A{2})*E(2);
A{3} = exp(P.A{3})*E(3);
C    = exp(P.C);
G    = exp(P.H)*diag(G);

% intrinsic connectivity and parameters
%--------------------------------------------------------------------------
Te   = T(1)/1000*exp(P.T(:,1));      % excitatory time constants
Ti   = T(2)/1000*exp(P.T(:,2));      % inhibitory time constants
Tk   = 512/1000;                     % slow potassium
He   = H(1)*exp(P.G);                % excitatory receptor density
Hi   = H(2);                         % inhibitory receptor density

% pre-synaptic inputs: s(V) with threshold adaptation
%--------------------------------------------------------------------------
R      = R.*exp(P.R);
x      = x';
X      = x;
X(1,:) = X(1,:) - X(13,:);
S      = 1./(1 + exp(-R(1)*(X - R(2)))) - 1./(1 + exp(R(1)*R(2)));
dSdx   = R(1)*exp(-R(1)*(max(X,-128) - R(2)))./(1 + exp(-R(1)*(X - R(2)))).^2;

% input
%==========================================================================
if isfield(M,'u')

    % endogenous input
    %----------------------------------------------------------------------
    U = u(:)*32;

else
    % exogenous input
    %----------------------------------------------------------------------
    U = C*u(:);
end

% State: f(x) and Jacobian dfdx
%==========================================================================

% NB: activity-dependent reduction in inhibitory effective time-constant
%--------------------------------------------------------------------------
Ti    = 4/1000 + Ti;

% intrinsic coupling
%--------------------------------------------------------------------------
for i = 1:n

    % synaptic dynamics - dfdx
    %----------------------------------------------------------------------
    ke    = -2/Te(i);
    ki    = -2/Ti(i);
    Ke    = -1/(Te(i)^2);
    Ki    = -1/(Ti(i)^2);
    dfdx  = [0   0   0   1   0   0   0   0   0   0   0   0   0
             0   0   0   0   1   0   0   0   0   0   0   0   0
             0   0   0   0   0   1   0   0   0   0   0   0   0
             Ke  0   0   ke  0   0   0   0   0   0   0   0   0
             0   Ke  0   0   ke  0   0   0   0   0   0   0   0
             0   0   Ki  0   0   ki  0   0   0   0   0   0   0
             0   0   0   0   0   0   0   1   0   0   0   0   0
             0   0   0   0   0   0  Ke   ke  0   0   0   0   0
             0   0   0   0   1   -1  0   0   0   0   0   0   0
             0   0   0   0   0   0   0   0   0   0   1   0   0
             0   0   0   0   0   0   0   0   0  Ki   ki  0   0
             0   0   0   0   0   0   0   1   0   0  -1   0   0
             0   0   0   0   0   0   0   0   0   0   0   0  -1/Tk];

    % intrinsic afferents - dfdS
    %----------------------------------------------------------------------
    Ke    = He(i)/Te(i);
    Ki    = Hi/Ti(i);

    j1    = Ke*G(i,1);
    j2    = Ke*G(i,2);
    j3    = Ke*G(i,3);
    j4    = Ki*G(i,4);
    j5    = Ki*G(i,5);
    dfdS  = [0   0   0   0   0   0   0   0   0   0   0   0   0
             0   0   0   0   0   0   0   0   0   0   0   0   0
             0   0   0   0   0   0   0   0   0   0   0   0   0
             0   0   0   0   0   0   0   0   j1  0   0   0   0
             j2  0   0   0   0   0   0   0   0   0   0   0   0
             0   0   0   0   0   0   0   0   0   0   0   j4  0
             0   0   0   0   0   0   0   0   0   0   0   0   0
             0   0   0   0   0   0   0   0   j3  0   0   0   0
             0   0   0   0   0   0   0   0   0   0   0   0   0
             0   0   0   0   0   0   0   0   0   0   0   0   0
             0   0   0   0   0   0   0   0   0   0   0  j5   0
             0   0   0   0   0   0   0   0   0   0   0   0   0
            4/Tk 0   0   0   0   0   0   0   0   0   0   0   0];

    % motion and Jacobian
    %----------------------------------------------------------------------
    dsdx       = diag(dSdx(:,i));
    dsdx(1,13) = -dsdx(1,1);
    dfdu       = sparse(4,1,Ke,s,1);

    F{i}       = dfdx*x(:,i) + dfdS*S(:,i) + dfdu*U(i);
    J{i,i}     = dfdx + dfdS*dsdx;

    % extrinsic afferents
    %----------------------------------------------------------------------
    for j = 1:n, if i ~= j

    k1    = Ke*(A{1}(i,j) + A{3}(i,j));
    k2    = Ke*(A{2}(i,j) + A{3}(i,j));
    k3    = Ki*(A{2}(i,j) + A{3}(i,j));
    dfdS  = [0   0   0   0   0   0   0   0   0   0   0   0   0
             0   0   0   0   0   0   0   0   0   0   0   0   0
             0   0   0   0   0   0   0   0   0   0   0   0   0
             0   0   0   0   0   0   0   0   k1  0   0   0   0
             0   0   0   0   0   0   0   0   k2  0   0   0   0
             0   0   0   0   0   0   0   0   0   0   0   0   0
             0   0   0   0   0   0   0   0   0   0   0   0   0
             0   0   0   0   0   0   0   0   k3  0   0   0   0
             0   0   0   0   0   0   0   0   0   0   0   0   0
             0   0   0   0   0   0   0   0   0   0   0   0   0
             0   0   0   0   0   0   0   0   0   0   0   0   0
             0   0   0   0   0   0   0   0   0   0   0   0   0
             0   0   0   0   0   0   0   0   0   0   0   0   0];

    % motion and Jacobian
    %----------------------------------------------------------------------
    F{i}   = F{i} + dfdS*S(:,j);
    J{i,j} = dfdS*diag(dSdx(:,j));

    end, end
end

% construct motion and Jacobian
%--------------------------------------------------------------------------
for i = 1:n
    k      = (1:n:s*n) + (i - 1);
    f(k,1) = F{i};
    for j  = 1:n
        l         = [1:n:s*n] + (j - 1);
        dfdx(k,l) = J{i,j};
    end
end

% extrinsic and intrinsic delays
%--------------------------------------------------------------------------
De = D(2).*exp(P.D)/1000;
Di = D(1).*exp(P.I)/1000;
De = (eye(n,n) - 1).*De;
Di = (eye(s,s) - 1)*Di;
De = kron(ones(s,s),De);
Di = kron(Di,eye(n,n));

D  = Di + De;

% Implement: dx(t)/dt = f(x(t + d)) = inv(1 - D.*dfdx)*f(x(t))
%--------------------------------------------------------------------------
D  = spm_inv(speye(n*s,n*s) - D.*dfdx);
f  = D*f;
J  = D*dfdx;

return

% Equations of motion
%==========================================================================

% Supragranular layer (inhibitory interneurons): depolarizing current
%--------------------------------------------------------------------------
f(:,7)  = x(:,8);
f(:,8)  = (He.*((A{2} + A{3})*S(:,9) + G(:,3).*S(:,9)) ...
           - 2*x(:,8) - x(:,7)./Te)./Te;

% Supragranular layer (inhibitory interneurons): hyperpolarizing current
%--------------------------------------------------------------------------
f(:,10) = x(:,11);
f(:,11) = (Hi*G(:,5).*S(:,12) ...
           - 2*x(:,11) - x(:,10)./Ti)./Ti;

% Granular layer (spiny stellate cells): depolarizing current
%--------------------------------------------------------------------------
f(:,1)  = x(:,4);
f(:,4)  = (He.*((A{1} + A{3})*S(:,9) + G(:,1).*S(:,9) + U) ...
           - 2*x(:,4) - x(:,1)./Te)./Te;

% Infra-granular layer (pyramidal cells): depolarizing current
%--------------------------------------------------------------------------
f(:,2)  = x(:,5);
f(:,5)  = (He.*((A{2} + A{3})*S(:,9) + G(:,2).*S(:,1)) ...
           - 2*x(:,5) - x(:,2)./Te)./Te;

% Infra-granular layer (pyramidal cells): hyperpolarizing current
%--------------------------------------------------------------------------
f(:,3)  = x(:,6);
f(:,6)  = (Hi*G(:,4).*S(:,12) ...
           - 2*x(:,6) - x(:,3)./Ti)./Ti;

% Surpa and Infra-granular layer (pyramidal cells): Voltage
%--------------------------------------------------------------------------
f(:,9)  = x(:,5) - x(:,6);
f(:,12) = x(:,8) - x(:,11);

% Granular layer (spiny stellate cells): hyperpolarizing current
%--------------------------------------------------------------------------
f(:,13) = (4*S(:,1) - x(:,13))./Tk;

% Jacobian for delays (evaluate numerically for simplicity)
%==========================================================================
dfdx = spm_diff('spm_fx_lfp',x,u,P,1,1);

author: spm


text:
function out = logdet(x)
a=chol(x);
out=sum(log(diag(a))*2);

end

author: european-central-bank

repo: BEAR-toolbox


text:
function plotps(x,Ts,N)
% Plot the power spectrum.
% plotps(x,[Ts,N])

if nargin < 2, Ts = 1; end
if nargin < 3, N  = 2^(ceil(log(length(x))/log(2))+2); end

Fs = 1/Ts;
X  = abs(fft(x,N)).^2;
%% X  = abs(2*pi*fft(x,N)/N*Ts).^2;

figure, plot(linspace(0,Fs/2,N/2+1),X(1:N/2+1));
set(gca,'FontSize',8,'XLim',[0 Fs/2],'YScale','log'); xlabel('Frequency'); ylabel('Amplitude');

%% figure, plot(linspace(0,Fs,N+1),[X,X(1)]);
%% set(gca,'FontSize',8,'XLim',[0 Fs],'YScale','log'); xlabel('Frequency'); ylabel('Amplitude');


text:
function SIIval = SII(E,N,Mtype)

% Implements  the ANSI S3.5-1997 standard:
% "Methods for calculation of the Speech Intelligibility Index".
%
%   INPUT:
%
%     'E' Speech Spectrum Level (Section 3.6 in the standard)
%         Level needs to be in dB SPL
%     'N' Noise Spectrum Level (Section 3.15 in the standard)
%         Level needs to be in dB SPL
%     'M' Speech material (needed to specify band-importance function)
%         A scalar having a value of either 1, 2, 3, 4, 5, 6. The Band-importance functions associated with each scalar are
%		                1:	various nonsense syllable tests where most English phonemes occur equally often (as specified in Table B.2)
%		                2:	CID-22 (as specified in Table B.2)
%		                3:	NU6 (as specified in Table B.2)
%		                3:	Diagnostic Rhyme test (as specified in Table B.2)
%		                5:	short passages of easy reading material (as specified in Table B.2)
%		                6:	SPIN (as specified in Table B.2)
%
%    OUTPUT:
%       Returns the SII value (0 to 1)
%
%     Note that only the one-third-octave band procedure is implemented.
%     Dimensions of 'E' and 'N' are 18x1, containing the speech and noise levels
%     at 1/3-octave frequencies (see line 78 for 1/3-octave center frequencies)
%
% EXAMPLE
% sp=[40 45 50 24 56 60 55 55 52 48 50 51 55 67 76 67 56 31];
% ns=[30 50 60 20 60 50 70 45 80 40 60 20 60 22 55 50 67 40];
% M= 5;
% sv = SII (sp,ns, M);
%
% (c)2012 Philipos C. Loizou

if length(E)~=18 | length(N)~=18
    error('The target and masker spectra vectors have incorrect dimension - needs to be 18.');
end
if Mtype>6 | Mtype<1
  error('Band-importance function type takes values between 1 and 6');
end

%================== DEFINE INPUT VARIABLES ==============================

G=zeros(1,18);  % insertion gains - needed if used in the context of amplification devices (hearing aids)
T=G;            % threshold levels in dB HL (for normal-hearing listeners, T=[0 0 ...0] )
VocalEffort = 'normal';  % or "raised", "loud" and "shout"

% Standard speech spectrum level for different vocal efforts (Table 3)
%
SpV=[32.41	33.81	35.29	30.77;
	34.48	33.92	37.76	36.65;
	34.75	38.98	41.55	42.5;
	33.98	38.57	43.78	46.51;
	34.59	39.11	43.3	47.4;
	34.27	40.15	44.85	49.24;
	32.06	38.78	45.55	51.21;
	28.3	36.37	44.05	51.44;
	25.01	33.86	42.16	51.31;
	23		31.89	40.53	49.63;
	20.15	28.58	37.7	47.65;
	17.32	25.32	34.39	44.32;
	13.18	22.35	30.98	40.8;
	11.55	20.15	28.21	38.13;
	9.33	16.78	25.41	34.41;
	5.31	11.47	18.35	28.24;
	2.59	7.67	13.87	23.45;
	1.13	5.07	11.39	20.72];

switch lower(VocalEffort)
	case 'normal', EV = SpV(:,1)';
	case 'raised', EV = SpV(:,2)';
	case 'loud',   EV = SpV(:,3)';
	case 'shout',  EV = SpV(:,4)';
	otherwise, error('Unknown level of vocal effort')
end;

% Define band center frequencies for 1/3rd octave procedure (Table 3)
f = [160 200 250 315 400 500 630 800 1000 1250 1600 2000, ...
     2500 3150 4000 5000 6300 8000];

% Define Internal Noise Spectrum Level (Table 3)
X = [0.6 -1.7 -3.9 -6.1 -8.2 -9.7 -10.8 -11.9 -12.5 -13.5 -15.4 -17.7, ...
	-21.2 -24.2 -25.9 -23.6 -15.8 -7.1];

% ----------------- start processing ----------------------
%
% Equivalent Speech Spectrum Level (5.1.3, Eq. 17)
E = E + G;

% Self-Speech Masking Spectrum (Sec 4.3.2.1, Eq. 5)
V = E - 24;

% 4.3.2.2
B = max(V,N+G);

% Calculate Equivalent Masking Spectrum Level (Sec 4.3.2.5, Eq. 9)
%
C = 0.6.*(B+10*log10(f)-6.353) - 80;  % slope parameter Ci (4.3.2.3 Eq. 7)
Z(1) = B(1);

for i = 2:18
	Z(i) = 10*log10(10.^(0.1*N(i)) + ...
	sum(10.^(0.1*(B(1:(i-1))+3.32.*C(1:(i-1)).*log10(0.89*f(i)./f(1:(i-1)))))));
end;

% Equivalent Internal Noise Spectrum Level (Sec 4.4 Eq. 10)
X = X + T;

% Compute disturbance Spectrum Level (4.5)
D = max(Z,X);

% Level Distortion Factor (Sec 4.6, Eq. 11)
%
L = 1 - (E - EV - 10)./160;
L = max(0,min(L,1));

% 4.7.1 Eq. 12
K = (E-D+15)/30;
K=max(0,min(K,1));

% Band Audibility Function (7.7.2 Eq. 13)
%
A = L.*K;

% Speech Intelligibility Index (4.8 Eq. 14)
%
SIIval = sum(BandImportance(Mtype).*A);

return;

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

function BIF = BandImportance(type)
%
% Band importance functions, taken from Table B.2:
% type  =
%		1:	Nonsense syllable tests where most English
%			phonemes occur equally often
%		2:	CID-22
%		3:	NU6 (monosyllables)
%		4:	Diagnostic Rhyme test (DRT)
%		5:	short passages of easy reading material
%		6:	SPIN (monosyllables)

if (nargin ~= 1) | (type>6) | (type<1)
	error('Incorrect argument to BandImportance');
end;

IFu   = [	0		0.0365	0.0168	0		0.0114	0
			0		0.0279	0.013	0.024	0.0153	0.0255
            0.0153  0.0405	0.0211	0.033	0.0179	0.0256
			0.0284	0.0500	0.0344	0.039	0.0558	0.036
			0.0363	0.0530	0.0517	0.0571	0.0898	0.0362
			0.0422	0.0518	0.0737	0.0691	0.0944	0.0514
			0.0509	0.0514	0.0658	0.0781	0.0709	0.0616
			0.0584	0.0575	0.0644	0.0751	0.066	0.077
			0.0667	0.0717	0.0664	0.0781	0.0628	0.0718
			0.0774	0.0873	0.0802	0.0811	0.0672	0.0718
			0.0893	0.0902	0.0987	0.0961	0.0747	0.1075
			0.1104	0.0938	0.1171	0.0901	0.0755	0.0921
			0.112	0.0928	0.0932	0.0781	0.082	0.1026
			0.0981	0.0678	0.0783	0.0691	0.0808	0.0922
			0.0867	0.0498	0.0562	0.048	0.0483	0.0719
			0.0728	0.0312	0.0337	0.033	0.0453	0.0461
			0.0551	0.0215	0.0177	0.027	0.0274	0.0306
			0		0.0253	0.0176	0.024	0.0145	0];

BIF = IFu(:,type)';

author: jtkim-kaist

repo: Speech-enhancement


text:
function ll = gpsimCandidateLogLikelihood(model)

% GPSIMCANDIDATELOGLIKELIHOOD Compute the log likelihood of a gene.
% FORMAT
% DESC computes the log likelihood of one gene given the existing genes
% in the GPSIM model. For use when ranking other prospective candidate genes.
% ARG model : the model containing the existing genes for which the log likelihood is computed.
% RETURN ll : the log likelihood of the new gene.
%
% SEEALSO : gpsimCreate, gpsimAddCandidate, gpsimCandidateLogLikelihood,
% gpsimCandidateObjective
%
% COPYRIGHT : Neil D. Lawrence, 2007

% SHEFFIELDML

dim = size(model.y, 1);

invK_ffm = model.candidate.invK*model.candidate.m;
K_ufinvK_ffm = model.candidate.K_uf*invK_ffm;

ll = -0.5*(dim*log(2*pi) + model.candidate.logDetK ...
           + model.candidate.logDetA ...
           - model.logDetK ...
           + sum(sum(invK_ffm.*model.candidate.m))...
           - sum(sum((model.candidate.Ainv*K_ufinvK_ffm).*K_ufinvK_ffm)));

% In case we need priors in.
if isfield(model.candidate, 'bprior'),
  ll = ll + kernPriorLogProb(model.candidate.kern);
  ll = ll + priorLogProb(model.candidate.bprior, model.B);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function visualize_adv(W, horiz, n_rows, n_cols, do_sort, borders)

if (nargin < 2)
    horiz = 0;
end

if (nargin < 5)
    do_sort = 0;
end

if (do_sort == 1)
    Wnorms = sum(W.^2,1);
    [B,IX] = sort(Wnorms,'descend');
    W = W(:,IX);
end

% how many pixels for borders?
if (nargin < 6)
    borders = 1;
end

fprintf('Visualizing patches\n');
[ndim,nunits]=size(W);
npix = floor(sqrt(ndim)+0.999);
npix2 = floor(sqrt(nunits)+0.999);
minW=min(W(:));
maxW=max(W(:));

if (nargin < 4)
    n_rows = npix2;
    n_cols = npix2;
end

%bigpic = -(minW+maxW)/2*ones(((npix+borders)*npix2+borders));
%bigpic = -(minW+maxW)/2*ones(((n_rows+borders)*npix2+borders), ((n_cols+borders)*npix2+borders));
bigpic = -(minW+maxW)/2*ones(((npix+borders)*n_rows+borders), ((npix+borders)*n_cols+borders));
%if (nunits/npix2<=npix2-1),
%    bigpic = bigpic(:,1:(npix+borders)*(npix2-1)+borders);
%end;
idx = 0;
for i=1:n_rows
    for j=1:n_cols
        idx = idx + 1;

        if idx > nunits
            break;
        end

        if (horiz)
            bigpic((i-1)*(npix+borders)+borders+1:(i-1)*(npix+borders)+borders+npix,...
                (j-1)*(npix+borders)+borders+1:(j-1)*(npix+borders)+borders+npix)...
                = reshape(W(:,idx),npix,npix)';
        else
            bigpic((i-1)*(npix+borders)+borders+1:(i-1)*(npix+borders)+borders+npix,...
                (j-1)*(npix+borders)+borders+1:(j-1)*(npix+borders)+borders+npix)...
                = reshape(W(:,idx),npix,npix);
        end
    end
end;
imagesc(bigpic);
colormap(gray);
axis off;
axis equal;
fprintf('done.\n');

author: kyunghyuncho

repo: deepmat


text:
function planC = importVarianFDF_3D(dirT,optFile)
% function planC = importVarianFDF_3D(dirT,optFile)
% Input: directory containing multi-frame FDF file and CERRoptions
% file or optS structure
% Output: CERR plan (planC)
% Example:
% dirT = 'C:\Projects\FDF reader\Test_data_Joel\3D-Large';
% optFile = 'CERROptions.m';
% planC = importVarianFDF_3D(dirT,optFile);
%
% APA, 05/15/2007
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

% Check the format of optFile.
if ~isstruct(optFile)
    optS = opts4Exe(optFile);
else
    optS = optFile;
end

planC = initializeCERR;
indexS = planC{end};
planC{indexS.CERROptions} = optS;
planC{indexS.indexS}     = indexS;

pathT=[dirT,'/'];

%%%% Import Transverse slices and interpolate

% 1> Read [X,Y,Z,f] from transverse slices

%%%%% Extract Data for Transverse Section
contents = dir(fullfile(pathT,'*.fdf'));
[header, data] = readVarianFDF_3D(pathT, contents(1).name);
for i=1:size(data,3)
    H{i}=header; DATA{i}=data(:,:,i);
    %%% Extract info from header
    fbr1=find(header=='{');
    fbr2=find(header=='}');
    MatrixDim(i,:)=str2num(header(fbr1(1)+1:fbr2(1)-1));
    span(i,:)=str2num(header(fbr1(4)+1:fbr2(4)-1));
    offset(i,:)=str2num(header(fbr1(5)+1:fbr2(5)-1));
    sliceCenter(i,:)=str2num(header(fbr1(8)+1:fbr2(8)-1));
    xC(i)=sliceCenter(i,1); yC(i)=sliceCenter(i,2);
    %zC(i)=sliceCenter(i,3);
    roi(i,:)=str2num(header(fbr1(9)+1:fbr2(9)-1));
    orientation(i,:)=str2num(header(fbr1(10)+1:fbr2(10)-1));
end
zC = linspace(0,span(1,3),MatrixDim(1,3));
sizeOfDimension1=MatrixDim(1,1);
sizeOfDimension2=MatrixDim(1,2);
sizeOfDimension1U=sizeOfDimension1; sizeOfDimension2U=sizeOfDimension2;
sizeDim1 = sizeOfDimension1-1;
sizeDim2 = sizeOfDimension2-1;
xSpan=span(1,1); ySpan=span(1,2);
xOffset=sliceCenter(1,1); yOffset=sliceCenter(1,2);
xOffsetU=xOffset; yOffsetU=yOffset;
grid2Units=xSpan/sizeDim2;
grid1Units=ySpan/sizeDim1;
voxelThickness=zC(2)-zC(1);

% xVals = xOffset - (sizeDim2*grid2Units)/2 : grid2Units : xOffset + (sizeDim2*grid2Units)/2;
% yVals = yOffset - (sizeDim1*grid1Units)/2 : grid1Units : yOffset + (sizeDim1*grid1Units)/2;
% zVals = zC;
% [Xt,Yt] = meshgrid(xVals,yVals);
% XXt=Xt(:); YYt=Yt(:);
% XTrans=[]; DD=[]; DDt=[];
% for i=1:length(zVals)
%     XTrans=[XTrans;
%         XXt, YYt, zVals(i)*XXt.^0];
%         DD1=DATA{i};
%     DD=[DD;
%         DD1(:)];
%     DDt=[DDt;
%         DD1(:)];
% end
%
% %% Reshape Data
% DDT3=reshape(DDt,[sizeOfDimension1,sizeOfDimension2,length(zVals)]);

DDT3 = data;

for i=1:length(zC)
    scanInfo(i).grid1Units=grid1Units;
    scanInfo(i).grid2Units=grid2Units;
    scanInfo(i).sizeOfDimension1=sizeOfDimension1U;
    scanInfo(i).sizeOfDimension2=sizeOfDimension2U;
    scanInfo(i).xOffset=xOffsetU;
    scanInfo(i).yOffset=yOffsetU;
    scanInfo(i).voxelThickness=voxelThickness;
    scanInfo(i).zValue=zC(i);
    scanInfo(i).CTOffset=1000;
    scanInfo(i).sliceThickness=[];
    scanInfo(i).imageType = 'MR Scan';
end
planC{3}(1).scanInfo = scanInfo;

minScanArray = min(DDT3(:));
maxScanArray = max(DDT3(:));
scanArray = (DDT3-minScanArray)/maxScanArray*4095;
planC{3}.scanArray = scanArray;
planC{3}.scanType = 'MRI';
planC{3}.scanUID = createUID('scan');
planC = setUniformizedData(planC,optS);

return

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function output = callfiltersdsp(model)

options = model.options.filtersd;

model = yalmip2nonlinearsolver(model);

if ~model.derivative_available
    disp('Derivate-free call to filterSD not yet implemented')
    error('Derivate-free call to filterSD not yet implemented')
end
if model.options.savedebug
    save filtersddebug model
end
showprogress('Calling filterSDSP',model.options.showprogress);

cu = [ repmat(0,length(model.bnonlinineq),1);
    repmat(0,length(model.bnonlineq),1);
    repmat(0,length(model.b),1);
    repmat(0,length(model.beq),1)];

cl = [ repmat(-inf,length(model.bnonlinineq),1);
    repmat(0,length(model.bnonlineq),1);
    repmat(-inf,length(model.b),1);
    repmat(0,length(model.beq),1)];

if isempty(cu)
    Flow = [];
    Fupp = [];
end

% These are needed to avoid recomputation due to ipopts double call to get
% f and df, and g and dg
global latest_x_f
global latest_x_g
global latest_df
global latest_f
global latest_G
global latest_g
global latest_xevaled
global latest_x_xevaled
latest_G = [];
latest_g = [];
latest_x_f = [];
latest_x_g = [];
latest_xevaled = [];
latest_x_xevaled = [];

funcs.objective = @(x)ipopt_callback_f(x,model);
funcs.gradient = @(x)ipopt_callback_df(x,model);
if ~isempty(cu)
    funcs.constraints = @(x)ipopt_callback_g(x,model);
    funcs.jacobian  = @(x)ipopt_callback_dg(x,model);
else
    funcs.constraints = [];
    funcs.jacobian = [];
end

lb  = model.lb(:)';
ub = model.ub(:)';

if ~isempty(cu)
   nljacstr = jacobiansparsityfromnonlinear(model);
else
   nljacstr = [];
end

if ~model.options.usex0
    model.x0 = (lb+ub)/2;
    model.x0(isinf(ub)) = lb(isinf(ub))+1;
    model.x0(isinf(lb)) = ub(isinf(lb))-1;
    model.x0(isinf(model.x0)) = 0;
end

options.display = model.options.verbose;
solvertime = tic;
[xout,fval,exitflag,stats,lambda] = filtersdsp(funcs.objective, funcs.gradient, model.x0, lb, ub, funcs.constraints, funcs.jacobian,nljacstr, cl, cu, options);
solvertime = toc(solvertime);

% Duals currently not supported
lambda = [];

x = RecoverNonlinearSolverSolution(model,xout);

switch exitflag
    case {0}
        problem = 0;
    case {1}
        problem = 2;
    case {2,3}
        problem = 1;
    case {5}
        problem = 3;
    case {4,9}
        problem = 4;
    case 105
        problem = 16;
    otherwise
        problem = -1;
end

% Internal format for duals
D_struc = [];

% Save all data sent to solver?
if model.options.savesolverinput
    solverinput.model = model;
else
    solverinput = [];
end

% Save all data from the solver?
if model.options.savesolveroutput
    solveroutput.x = xout;
    solveroutput.exitflag = exitflag;
    solveroutput.stats = stats;
else
    solveroutput = [];
end

% Standard interface
output = createoutput(x,D_struc,[],problem,'filterSD',solverinput,solveroutput,solvertime);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [speciespKas,includedMets] = assignpKasToSpecies(mets,neutralMolfileDir,msDistrDir,pKaDir)

% Use each metabolite's InChI and ChemAxon's cxcalc (calculator plugin) to
% get pKas(I=0), charges and numbers of protons for all microspecies
% present at pH 5-9.
%
% [speciespKas,includedMets] = assignpKasToSpecies(mets,neutralMolfileDir,msDistrDir,pKaDir)
%
% INPUTS
% mets                  Cell array of metabolite ID. Current implementation
%                       assumes ID are formatted as BiGG_ID[Compartment];
%                       e.g., atp[c] for ATP in cytosol.
% neutralMolfileDir     Path to directory containing molfiles for neutral
%                       metabolites. Molfiles must be named with the
%                       metabolite ID in mets; e.g., atp.mol.
% msDistrDir            Path to directory containing metabolite species
%                       distributions at pH 5, 6, 7, 8 and 9. Species
%                       distributions are created with the ChemAxon
%                       calculator plugin. Names of SD files containing
%                       species distributions should be formatted as
%                       metID_msdistr_pHA.sdf; e.g., atp_msdistr_pH5.sdf.
% pKaDir                Path to directory containing pKa estimates for
%                       metabolites. pKas are estimated with the ChemAxon
%                       calculator plugin. Names of text files containing
%                       pKa estimates should be formatted as
%                       metID_pkas.txt; e.g., atp_pkas.txt.
%
% OUTPUTS
% speciespKas       Structure containing the following fields for each
%                   metabolite:
% .inchis           n x 1 cell array with a species-specific InChI string
%                   for each species. inchis(1) = InChI for species 1 etc.
%                   NOTE: Species 1 is the species with the fewest hydrogen
%                   atoms etc.
% .formulas         n x 1 cell array with species-specific chemical
%                   formulas.
% .zs               n x 1 vector containing the charge on each species.
% .nHs              n x 1 vector containing the number of hydrogen atoms in
%                   each species.
% .majorMSpH7       n x 1 boolean vector with a logical 1 in the row for
%                   the major microspecies at pH 7 according to ChemAxon.
% .abundanceAtpH7   n x 1 vector with species percentage abundance at pH 7.
% .pKas             n x n matrix where element i,j is the pKa for the
%                   acid-base equilibrium between species i and j.
% includedMets      Cell array of metabolite ID for those metabolites that
%                   were included in speciespKas. Should include all
%                   metabolites that had molfiles in neutralMolfileDir.

phs = 5:9; % Physiological pH range

dirData = dir(neutralMolfileDir);
molfileList = cat(2,dirData.name);
molfileList = regexp(molfileList,'.mol','split');
molfileList = regexprep(molfileList,'\.','');
molfileList = molfileList(~cellfun('isempty',molfileList));

includedMets = {};

model.mets = mets;
nMets = length(model.mets);
done = cell(nMets,1);

for k = 1:nMets
    disp(['Starting on metabolite ' num2str(k) '. Only ' num2str(nMets-k) ' metabolites to go.'])

    if ~any(strcmp(model.mets{k}(1:end-3),done))
        done{k} = model.mets{k}(1:end-3);

        if any(strcmp(molfileList,model.mets{k}(1:(end-3))))
            includedMets = [includedMets; model.mets{k}(1:end-3)];

            % Get data from sdf files with microspecies distributions
            % for all species with molfractions above threshold
            totalCharges = {};
            chargedAtoms = {};
            atomicCharges = {};
            currentPhs = {};
            percentages = {};
            mmsPh7Bool = {};

            fid5 = fopen('abundantSpecies.sdf', 'w+'); % Create separate sdf file for species with molfractions above threshold
            fclose(fid5);

            for n = 1:length(phs)
                ph = phs(n);

                newPh = 1; % True while processing first microspecies at each pH

                %     Open microspecies distribution at ph
                fid2 = fopen([msDistrDir filesep model.mets{k}(1:(end-3)) '_msdistr_ph' num2str(ph) '.sdf'], 'r'); % Open sdf file with microspecies distribution for reading

                charged = 0; % Neutral until proven charged

                %     Parse speciesDistribution.sdf for all microspecies that are present at ph
                while 1
                    tline = fgetl(fid2); % Get next line of speciesPhDistribution.sdf

                    if ~ischar(tline) % Exit while loop if the last line of speciesPhDistribution.sdf is reached, continue at next pH
                        fclose(fid2);
                        delete currentSpecies.sdf
                        pause(eps);
                        break;
                    end

                    fid3 = fopen('currentSpecies.sdf', 'a'); % Create separate sdf file for species with molfractions above threshold
                    fprintf(fid3, '%s\r\n', tline); % Copy line to abundantSpecies.sdf
                    fclose(fid3);

                    if ~isempty(strfind(tline, 'M  CHG'))
                        charged = 1;
                        chargeLine = tline;
                    end

                    if ~isempty(strfind(tline, '>  <DISTR')) % Get the percentage of current microspecies at ph
                        percentLine = fgetl(fid2);
                        fid3 = fopen('currentSpecies.sdf', 'a'); % Create separate sdf file for species with molfractions above threshold
                        fprintf(fid3, '%s\r\n', percentLine); % Copy line to abundantSpecies.sdf
                        fclose(fid3);

                        percentage = str2double(regexprep(percentLine, ',', '\.'));

                        if percentage == 0 % Exit while loop at first microspecies that is not present at current pH. Continue at next pH.
                            fclose(fid2);
                            delete currentSpecies.sdf
                            pause(eps);
                            break;
                        else
                            currentPhs = [currentPhs; {ph}];
                            percentages = [percentages; {percentage}];

                            % Get numbers of atoms carrying charge, charge on each atom
                            % and total charge
                            currentChargedAtoms = [];
                            currentAtomicCharges = [];

                            if charged
                                splitLine = regexp(strtrim(chargeLine(11:end)), '\s+', 'split');
                                AtomIdx = 1:2:(length(splitLine)-1);
                                chargeIdx = 2:2:length(splitLine);

                                for m = 1:length(AtomIdx)
                                    currentChargedAtoms = [currentChargedAtoms, str2double(splitLine{AtomIdx(m)})];
                                end

                                for m = 1:length(chargeIdx)
                                    currentAtomicCharges = [currentAtomicCharges, str2double(splitLine{chargeIdx(m)})];
                                end

                                currentTotalCharge = sum(currentAtomicCharges);

                            else
                                currentTotalCharge = 0;

                            end

                            chargedAtoms = [chargedAtoms; {currentChargedAtoms}];
                            atomicCharges = [atomicCharges; {currentAtomicCharges}];
                            totalCharges = [totalCharges; {currentTotalCharge}];
                        end
                    end

                    if ~isempty(strfind(tline, '$$$$')) % End of sdf for current species
                        system('type currentSpecies.sdf >> abundantSpecies.sdf');
                        delete currentSpecies.sdf
                        pause(eps);
                        if ph == 7 && newPh
                            mmsPh7Bool = [mmsPh7Bool; {true}]; % True for major microspecies at pH 7, the one WebGCM computes Gibbs energy of formation for
                        else
                            mmsPh7Bool = [mmsPh7Bool; {false}];
                        end

                        newPh = 0; % reset
                        charged = 0; % Reset
                    end
                end
            end

            % Get additional species data
            % Get formulas
            [~, formulas] = system('cxcalc formula abundantSpecies.sdf');
            formulas = regexp(formulas,'\s','split');
            formulas = formulas(4:2:2*length(mmsPh7Bool)+3)';

            % Get number of hydrogen atoms
            numberHs = cell(size(formulas));
            for spnr = 1:length(formulas)
                numberHs{spnr} = numAtomsOfElementInFormula(formulas{spnr},'H');
            end

            % Get inchis for all species in abundantSpecies.sdf
            [~, inchis] = system('babel abundantSpecies.sdf -oinchi -xFwT/nostereo');
            inchis = regexp(inchis, '\n', 'split');
            inchiIdx = regexp(inchis,'InChI=');
            inchis = inchis(~cellfun('isempty',inchiIdx))';
            inchiIdx = inchiIdx(~cellfun('isempty',inchiIdx))';

            if ~isempty(setdiff(1:length(inchis),strmatch('InChI=',inchis)))
                errlins = setdiff(1:length(inchis),strmatch('InChI=',inchis));
                for errlini = 1:length(errlins)
                    inchis{errlins(errlini)} = inchis{errlins(errlini)}(inchiIdx{errlins(errlini)}:end);
                end
            end

            % Create dataArray with all species data (non-unique)
            dataArray = [inchis formulas numberHs totalCharges chargedAtoms atomicCharges currentPhs percentages mmsPh7Bool];

            % Sum up percentages for different stereoisomers of the same
            % species at the same pH
            stiCorrectedDataArray = [];

            for stph = 5:9
                currentpHBool = cat(1,dataArray{:,7}) == stph;
                dataArrayCurrentpH = dataArray(currentpHBool,:);
                uniqueInchisCurrentpH = unique(inchis(currentpHBool));
                for stinr = 1:length(uniqueInchisCurrentpH)
                    stiCorrectedDataArray = [stiCorrectedDataArray; dataArrayCurrentpH(find(ismember(dataArrayCurrentpH(:,1),uniqueInchisCurrentpH{stinr}),1),:)];
                    stiCorrectedDataArray{end,8} = sum(cat(1,dataArrayCurrentpH{ismember(dataArrayCurrentpH(:,1),uniqueInchisCurrentpH{stinr}),8}));
                end
            end

            dataArray = stiCorrectedDataArray;

            % Discard species that are present below some logical
            % threshold percentage
            thresholdPercentage = 1;
            belowThresholdBool = cat(1,dataArray{:,8}) < thresholdPercentage;
            dataArray = dataArray(~belowThresholdBool,:);

            abundanceAtpH7 = zeros(size(dataArray,1),1);
            abundanceAtpH7(cat(1,dataArray{:,7}) == 7) = cat(1,dataArray{cat(1,dataArray{:,7})==7,8});
            dataArray = [dataArray, num2cell(abundanceAtpH7)];

            % Get inchis for all species in abundantSpecies.sdf
            %                 [~, uniqueInchis] = system('babel abundantSpecies.sdf -oinchi -xFwuT/nostereo');
            %                 uniqueInchis = regexp(uniqueInchis, '\n', 'split');
            %                 uniqueInchis = uniqueInchis(strmatch('InChI=',uniqueInchis))';
            inchis = dataArray(:,1);
            uniqueInchis = unique(inchis);

            uniqueDataArray = {};
            for inr = 1:length(uniqueInchis)
                sameSpeciesBool = ismember(dataArray(:,1),uniqueInchis{inr});
                tmpDataArray = dataArray(sameSpeciesBool,:);
                [~, speciesIdxInTmpDataArray] = max(cat(1,tmpDataArray{:,8}));
                uniqueDataArray = [uniqueDataArray; tmpDataArray(speciesIdxInTmpDataArray,:)];
                if any(cat(1,dataArray{sameSpeciesBool,9}))
                    uniqueDataArray{end,9} = true;
                end
                uniqueDataArray{end,10} = max(cat(1,dataArray{sameSpeciesBool,10}));
            end

            dataArray = uniqueDataArray;

            % Get charge and number of H atoms for fully protonated
            % species for later corrections
            nmoldir = fopen([neutralMolfileDir filesep model.mets{k}(1:(end-3)), '.mol'], 'r');
            nmolcharged = false;

            while 1
                tline = fgetl(nmoldir);

                if ~ischar(tline)
                    break;
                end

                if ~isempty(strfind(tline, 'M  CHG'))
                    nmolcharged = true;
                    nmolAcharges = regexp(tline(12:end),'\s+','split');
                    nmolAcharges = nmolAcharges(~cellfun('isempty',nmolAcharges));
                    nmolAcharges = nmolAcharges(2:2:end);
                    acnum = [];
                    for acIdx = 1:length(nmolAcharges)
                        acnum = [acnum, str2double(nmolAcharges(acIdx))];
                    end
                    nmolcharge = sum(acnum);
                end
            end

            fclose(nmoldir);

            if ~nmolcharged
                nmolcharge = 0;
            end

            [~, nmolformula] = system(['cxcalc formula ' neutralMolfileDir filesep model.mets{k}(1:(end-3)), '.mol']);
            nmolformula = regexp(nmolformula,'\s','split');
            nmolformula = nmolformula{4};
            nmolnH = numAtomsOfElementInFormula(nmolformula,'H');

            dataArray(:,3) = num2cell(nmolnH + (cat(1,dataArray{:,4}) - nmolcharge));
            dataArray = sortrows(dataArray,[3, 4]); % Sort by number of hydrogen atoms, then charge, in ascending order so row 1 contains data for species 1.

            % Boolean vector indicating major microspecies at pH 7
            mmsBool = logical(cat(1,dataArray{:,9}));

            % % Figure out which groups are in equilibrium at the specified pH by
            % comparing numbering of atoms carrying charge in each microspecies.
            nSpecies = size(dataArray,1);
            metIdx = length(includedMets);

            speciespKas(metIdx).abbreviation = includedMets{metIdx};
            speciespKas(metIdx).inchis = dataArray(:,1);
            speciespKas(metIdx).formulas = dataArray(:,2);
            speciespKas(metIdx).zs = cat(1,dataArray{:,4});
            speciespKas(metIdx).nHs = cat(1,dataArray{:,3});
            speciespKas(metIdx).majorMSpH7 = mmsBool;
            speciespKas(metIdx).abundanceAtpH7 = cat(1,dataArray{:,10});

            if nSpecies >= 2 % If there is only one species there is no need to compute pKas

                % Get species pKas
                maxChargedAtom = max(unique(cat(2,dataArray{:,5}))); % Find highest atomic number carrying charge
                chargeMatrix = zeros (nSpecies,maxChargedAtom); % Each column is an atom, each row a species. Each element is the charge on an atom in a species.
                atomNumbers = {};

                for n = 1:nSpecies
                    if ~isempty(dataArray{n,5}) % If empty the species is neutral and no atoms carry charge
                        chargeMatrix(n,dataArray{n,5}) = dataArray{n,6}; % Populate chargeMatrix
                    end
                end

                speciesEquilibriumBool = false(nSpecies,nSpecies);
                for spRows = 1:nSpecies
                    for spCols = 1:nSpecies
                        if length(find(chargeMatrix(spRows,:) ~= chargeMatrix(spCols,:))) == 1 && abs(abs(dataArray{spRows,4}) - abs(dataArray{spCols,4})) == 1
                            speciesEquilibriumBool(spRows,spCols) = true;
                        end
                    end
                end

                speciesEquilibriumAtomNrs = zeros(size(speciesEquilibriumBool));
                for n = 1:(nSpecies-1)
                    theseAtomNumbers = [];
                    oneup = find(speciesEquilibriumBool(n,:) & (1:nSpecies > n)); % Find row indices of all species with one more H atom

                    if ~isempty(oneup) % Not for species with the highest number of H atoms
                        for l = 1:length(oneup)
                            theseAtomNumbers = [theseAtomNumbers, find(chargeMatrix(n,:) ~= chargeMatrix(oneup(l),:))]; % Find number of atom being protonated/deprotonated in acid-base equilibrium between species n and species n+1
                            speciesEquilibriumAtomNrs(n,oneup(l)) = find(chargeMatrix(n,:) ~= chargeMatrix(oneup(l),:));
                        end

                        atomNumbers = [atomNumbers; {theseAtomNumbers}];
                    end
                end
                speciesEquilibriumAtomNrs = speciesEquilibriumAtomNrs + speciesEquilibriumAtomNrs';

                fid4 = fopen([pKaDir filesep model.mets{k}(1:(end-3)) '_pkas.txt'], 'r'); % Open txt file with pkas for reading
                headerline = fgetl(fid4);
                pkaline = fgetl(fid4);
                fclose(fid4);

                %     Construct two column cell array with pkas of all functional groups
                %     and numbers of atoms associated with each pka
                splitline = regexp(pkaline, '\s', 'split');
                pkaCells = splitline(2:end-1);
                pkaCells = pkaCells(~cellfun('isempty',pkaCells));
                pkaCells = regexprep(pkaCells, ',', '\.');

                splitAtomNumbers = regexp(splitline{end}, ',', 'split');
                pkaCells = [pkaCells;splitAtomNumbers]';
                pkaMat = zeros(size(pkaCells));

                for n = 1:size(pkaMat,1);
                    pkaMat(n,1) = str2double(pkaCells{n,1});
                    pkaMat(n,2) = str2double(pkaCells{n,2});

                end

                %     Correlate numbers of atoms associated with each pka with numbers of
                %     atoms in atomNumbers to get pka for equilibrium between each pair of
                %     microspecies
                speciesEquilibriumPkas = zeros(size(speciesEquilibriumBool));

                for n = 1:length(atomNumbers)
                    for m = 1:length(atomNumbers{n})
                        pkaIdx = find(atomNumbers{n}(m) == pkaMat(:,2));
                        if length(pkaIdx) == 1
                            speciesEquilibriumPkas(speciesEquilibriumAtomNrs == pkaMat(pkaIdx,2)) = pkaMat(pkaIdx,1);

                        elseif length(pkaIdx) > 1 % Many atoms will have two pKas; one acidic (for when it donates a proton) and one basic (for when it accepts a proton)
                            [~, minIdx] = min(abs(dataArray{n,7} - pkaMat(pkaIdx,1))); % Assume the appropriate pKa is the one closer to the pH where species is present at its maximum
                            speciesEquilibriumPkas(speciesEquilibriumAtomNrs == pkaMat(pkaIdx(minIdx),2)) = pkaMat(pkaIdx(minIdx),1);
                        end

                    end

                end

                speciespKas(metIdx).pKas = speciesEquilibriumPkas;

            else
                speciespKas(metIdx).abbreviation = includedMets{metIdx};
                speciespKas(metIdx).inchis = dataArray(:,1);
                speciespKas(metIdx).formulas = dataArray(:,2);
                speciespKas(metIdx).zs = cat(1,dataArray{:,4});
                speciespKas(metIdx).nHs = cat(1,dataArray{:,3});
                speciespKas(metIdx).majorMSpH7 = mmsBool;
                speciespKas(metIdx).abundanceAtpH7 = 100;
                speciespKas(metIdx).pKas = [];
            end

        end

    end

end

speciespKas = speciespKas';

author: opencobra

repo: cobratoolbox


text:
function arcsinh_values_test ( )

%*****************************************************************************80
%
%% ARCSINH_VALUES_TEST demonstrates the use of ARCSINH_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 June 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ARCSINH_VALUES_TEST:\n' );
  fprintf ( 1, '  ARCSINH_VALUES stores values of \n' );
  fprintf ( 1, '  the hyperbolic arc sine function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X           FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = arcsinh_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function plot_components_3D_GUI(Y,A,C,b,f,Cn,options)

% GUI for plotting components for the case of 3D volumetric imaging

memmaped = isobject(Y);
defoptions = CNMFSetParms;
if nargin < 7 || isempty(options); options = []; end
if ~isfield(options,'d1') || isempty(options.d1); d1 = input('What is the total number of rows? \n'); else d1 = options.d1; end          % # of rows
if ~isfield(options,'d2') || isempty(options.d2); d2 = input('What is the total number of columns? \n'); else d2 = options.d2; end          % # of columns
if ~isfield(options,'d3') || isempty(options.d2); d3 = input('What is the total number of z-stacks? \n'); else d3 = options.d3; end          % # of columns
if ~isfield(options,'plot_df') || isempty(options.plot_df); options.df = defoptions.plot_df; end
plot_df = options.plot_df;
if ~isfield(options,'full_A') || isempty(options.full_A); full_A = defoptions.full_A; else full_A = options.full_A; end

T = size(C,2);
if ~ismatrix(Y)
    Y = reshape(Y,d1*d2*d3,T);
end
if nnz(A)/numel(A) > 0.3; A = full(A); end
center = com(A,d1,d2,d3);
if size(center,2) == 2
    center(:,3) = 1;
end
center = round(center);
b = double(b);
C = double(C);
f = double(f);
nA = full(sqrt(sum(A.^2))');
[K,~] = size(C);
A = A/spdiags(nA,0,K,K);    % normalize spatial components to unit energy
C = bsxfun(@times,C,nA(:));

nb = size(f,1);     % number of background components

step = 5e3;
if memmaped
    AY = zeros(K,T);
    for j = 1:step:d
        AY = AY + A(j:min(j+step-1,d),:)'*double(Y.Yr(j:min(j+step-1,d),:));
    end
else
    if issparse(A) && isa(Y,'single')
        if full_A
            AY = full(A)'*Y;
        else
            AY = A'*double(Y);
        end
    else
        AY = A'*Y;
    end
end
Y_r = (AY- (A'*A)*C - full(A'*double(b))*f) + C;

if plot_df
    [~,Df] = extract_DF_F(Y,[A,double(b)],[C;f],[],options);
else
    Df = ones(size(A,2)+1,1);
end

fig = figure('Visible','off');
set(gcf,'Position',2*[300,300,960,480]);
set(gcf,'PaperPosition',2*[300,300,960,480]);

thr = 0.95;
minC = max(squeeze(min(min(Cn,[],1),[],2)),0);
maxC = squeeze(max(max(Cn,[],1),[],2));

% Create a figure and axes
% ax = axes('Units','DF/F');
% Create slider

sld = uicontrol('Style', 'slider',...
    'Min',1,'Max',K+nb,'Value',1,'SliderStep',[1/(K+nb-1) 1],...
    'Position', [150 20 800 20],...
    'Callback', @surfzlim);

% Add a text uicontrol to label the slider.
txt = uicontrol('Style','text',...
    'Position',[400 45 120 20],...
    'String','Component');

% Make figure visble after adding all components
fig.Visible = 'on';
plot_component(1)

% This code uses dot notation to set properties.
% Dot notation runs in R2014b and later.
% For R2014a and earlier: set(f,'Visible','on');

    function surfzlim(source,callbackdata)
        i = source.Value;
        plot_component(round(i))
        % For R2014a and earlier:
        % i = get(source,'Value');
    end

    function plot_component(i)
        if i <= K
            %subplot(3,2,5);
             subplot(3,3,7);
            cla
            Atemp = reshape(full(A(:,i)),d1,d2,d3);
            data = smooth3(Atemp);
            patch(isocaps(data,0.25*max(data(:))),...
               'FaceColor','interp','EdgeColor','none');
            p1 = patch(isosurface(data,0.25*max(data(:))),...
               'FaceColor','blue','EdgeColor','none');
            isonormals(data,p1)
            view(3);
            axis vis3d tight
            camlight left;
            %colormap jet
            lighting gouraud
            title(sprintf('Component %i',i),'fontsize',16,'fontweight','bold'); drawnow;
        end
        if i <= K
            %subplot(3,2,5);
             subplot(3,3,[2,5,8]);
            %cla
            Atemp = reshape(full(A(:,i)),d1,d2,d3);
            data = smooth3(Atemp);
            patch(isocaps(data,0.25*max(data(:))),...
               'FaceColor','interp','EdgeColor','none');
            p1 = patch(isosurface(data,0.25*max(data(:))),...
               'FaceColor','blue','EdgeColor','none');
            isonormals(data,p1)
            view(3);
            axis vis3d %tight
            set(gca,'XLim',[1,d1],'YLim',[1,d2],'ZLim',[1,d3]);
            camlight left;
            %colormap jet
            lighting gouraud
            title(sprintf('Components 1 - %i',i),'fontsize',16,'fontweight','bold'); drawnow;
        end
        %subplot(3,2,[1,3]);
        subplot(3,3,[1,4]);
        if i <= K
            cla
            %,[minC(center(i,3)),maxC(center(i,3))]
            imagesc(Cn(:,:,center(i,3)),[minC(center(i,3)),maxC(center(i,3))]); axis equal; axis tight; axis off; hold on; colorbar;
            A_temp = reshape(full(A(:,i)),d1,d2,d3);
            A_temp = A_temp(:,:,center(i,3));
            A_temp = medfilt2(A_temp,[3,3]);
            A_temp = A_temp(:);
            [temp,ind] = sort(A_temp(:).^2,'ascend');
            temp =  cumsum(temp);
            ff = find(temp > (1-thr)*temp(end),1,'first');
            if ~isempty(ff)
                [~,ww] = contour(reshape(A_temp,d1,d2),[0,0]+A_temp(ind(ff)),'LineColor','k');
                ww.LineWidth = 2;
            end
            title(sprintf('Component %i, z-plane %i',i,center(i,3)),'fontsize',16,'fontweight','bold'); drawnow; %pause;
        else
            cla
            b_temp = reshape(b(:,i-K),d1,d2,d3);
            imagesc(b_temp(:,:,ceil(d3/2))); axis equal; axis tight;
            title(['Background component, z-plane',num2str(ceil(d3/2))],'fontsize',16,'fontweight','bold'); drawnow;
        end
        %subplot(3,2,[2,4,6]);
        subplot(3,3,[3,6,9]);
        if i <= K
            plot(1:T,Y_r(i,:)/Df(i),'linewidth',2); hold all; plot(1:T,C(i,:)/Df(i),'linewidth',2);
            if plot_df
                title(sprintf('Component %i (calcium DF/F value)',i),'fontsize',16,'fontweight','bold');
            else
                title(sprintf('Component %i (calcium raw value)',i),'fontsize',16,'fontweight','bold');
            end
            leg = legend('Raw trace (filtered)','Inferred');
            set(leg,'FontSize',14,'FontWeight','bold');
            title(sprintf('Component %i',i),'fontsize',16,'fontweight','bold'); drawnow; %pause;
            drawnow;
            hold off;
        else
            plot(1:T,f(i-K,:)); title('Background activity','fontsize',16,'fontweight','bold');
            drawnow;
        end
    end
end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function [cplx_str,cplx_inc,cplx_lib,cplx_libname,cplx_ver] = opti_FindCplex()
%Finds IBM ILOG CPLEX Libraries and returns path and version

%Known CPLEX path locations (Modify to suit your system by adding to cell arrays, or create a new structure for other versions)
CPLX125.x64 = {'C:\Program Files\IBM\ILOG\CPLEX_Studio125\cplex\'};
CPLX125.x86 = {'C:\Program Files (x86)\IBM\ILOG\CPLEX_Studio125\cplex\'};
CPLX125.ver = '12.5';
CPLX125.lib = 'cplex125';

CPLX124.x64 = {'C:\Program Files\IBM\ILOG\CPLEX_Studio_Academic124\cplex\'};
CPLX124.x86 = {'C:\Program Files (x86)\IBM\ILOG\CPLEX_Studio_Academic124\cplex\'};
CPLX124.ver = '12.4';
CPLX124.lib = 'cplex124';

%Add any new structures to below
CPLXLIB = {CPLX125,CPLX124};

%DO NOT MODIFY BELOW HERE
for i = 1:length(CPLXLIB)
    %Check for CPLEX
    [ok,cplx_str,cplx_inc,cplx_lib,cplx_libname,cplx_ver] = checkCPLXVer(CPLXLIB{i});
    if(ok), return; end
end

%Still no Luck - user must locate
error('Could not find the IBM ILOG CPLEX Location on your computer. Please modify this file to locate it');

function [ok,cplx_str,cplx_inc,cplx_lib,cplx_libname,cplx_ver] = checkCPLXVer(cplxstr)
%Local check function
ok = false;
%Defaults
cplx_str = []; cplx_inc = []; cplx_lib = []; cplx_libname = []; cplx_ver = [];
%Find local config paths
switch(computer)
    case 'PCWIN'
        if(~isfield(cplxstr,'x86') || isempty(cplxstr.x86))
            error('No x86 CPLEX directories specified to check');
        end
        paths = cplxstr.x86;
    case 'PCWIN64'
        if(~isfield(cplxstr,'x64') || isempty(cplxstr.x86))
            error('No x64 CPLEX directories specified to check');
        end
        paths = cplxstr.x64;
    otherwise
        error('This function is only setup for Windows PCs');
end
for i = 1:length(paths)
    if(exist(paths{i},'dir'))
        dir_cplx = paths{i};
        %Get Include Directory
        cplx_inc = [dir_cplx 'include\ilcplex'];
        if(~exist(cplx_inc,'dir'))
            error('Could not find Cplex include directory. Checked:\n %s',cplx_inc);
        end
        %Get Library Directory
        switch(computer)
            case 'PCWIN'
                cplx_lib = [dir_cplx 'lib\x86_windows_vs2010\stat_mda\'];
            case 'PCWIN64'
                cplx_lib = [dir_cplx 'lib\x64_windows_vs2010\stat_mda\'];
        end
        if(~exist(cplx_lib,'dir'))
            error('Could not find Cplex lib directory. Checked:\n %s\n',cplx_lib);
        end
        %Build complete linker line
        cplx_str = [' -I"' cplx_inc '" -L"' cplx_lib '" -l' cplxstr.lib ' '];
        cplx_libname = cplxstr.lib;
        cplx_ver = cplxstr.ver;
        ok = true;
        return;
    end
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [  ] = NiceSave(figname,figfolder,baseName,varargin)
%NiceSave(figname,figfolder,baseName,figtype) formats the figure for best viewing
%and saves it as a .pdf in figfolder with name recname_figname.pdf
%
%INPUTS
%   figname     string with title of figure
%   figfolder   folder to save figure to
%   baseName    name of the recording ([] for none)
%   'figtype'     'pdf','jpeg',etc
%   'includeDate' default: false. include date in filename
%
%DLevenstein Fall 2016
%TODO: add option to not overwrite old figure, but instead save as
%'figname02' etc.
%%
p = inputParser;
addParameter(p,'figtype','pdf')
addParameter(p,'includeDate',false)
parse(p,varargin{:})
figtype = p.Results.figtype;
includeDate = p.Results.includeDate;

%%

if ~exist(figfolder,'dir')
    mkdir(figfolder)
end

if includeDate
   figname = [figname,'_',date];
end
%set(gcf,'TickDir','out')

%if strcmp(figtype,'jpg')
set(gcf,'PaperOrientation','landscape');
%end
set(gcf,'PaperUnits','normalized');
set(gcf,'PaperPosition',[0 0 1 1]);
%orient portrait
saveas(gcf,[figfolder,'/',baseName,'_',figname,'.',figtype],figtype) ;

end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function Tbig = divide_by_pot(Tbig, Tsmall)
% DIVIDE_BY_POT Tbig /= Tsmall
% Tbig = divide_by_pot(Tbig, Tsmall)
%
% Tsmall's domain must be a subset of Tbig's domain.

smallp = extend_domain_table(Tsmall.p, Tsmall.domain, Tsmall.sizes, Tbig.domain, Tbig.sizes);
smallp = smallp + (smallp==0);
Tbig.p = Tbig.p ./ smallp;

smallu = extend_domain_table(Tsmall.u, Tsmall.domain, Tsmall.sizes, Tbig.domain, Tbig.sizes);
Tbig.u = Tbig.u - smallu;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt