text:
function [ po, pc, pe ] = lagrange_partial2 ( d, n, r, nd, xd )

%*****************************************************************************80
%
%% LAGRANGE_PARTIAL2: Partial Lagrange polynomial basis from data.
%
%  Discussion:
%
%    This function represents algorithm 4.1 in the reference,
%    modified for the case where the number of data points is less
%    than the dimension of the desired polynomial space,
%    with the further modification that a form of "pivoting" is used
%    to select the next polynomial as the one with maximum absolute
%    value at the current node.
%
%    This function is given XD, a set of ND distinct data points in a
%    D dimensional space, and returns information defining a set of
%    ND Lagrange polynomials L(i)(X) with the property that:
%
%      L(i)(XD(j)) = delta(i,j)
%
%    This function is used in cases where ND, the number of data points,
%    is less than or equal to R, the dimension of the space of polynomials
%    in D dimensions and total degree N or less, that is:
%
%      ND <= R = Choose ( N + D, N )
%
%    There will be ND polynomials returned.  Each polynomial can have
%    as many as R coefficients.
%
%    Each polynomial is given as a vector, with each entry corresponding
%    to a nonzero coefficient.  In particular, for polynomial L(i)(X):
%
%      PO(i) is the order, that is, the number of nonzero coefficients;
%      PC(i,j), for 1 <= j <= PO(i), is the coefficient of the J-th term.
%      PE(i,j), for 1 <= j <= PO(i), encodes the exponents of the J-th term.
%
%    The exponent codes are a compact way of recording the exponent vector
%    associated with each monomial.  If PE(i,j) = k, then the corresponding
%    vector of D exponents can be determined by:
%
%      E = mono_unrank_grlex ( D, k );
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 January 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Tomas Sauer, Yuan Xu,
%    On multivariate Lagrange interpolation,
%    Mathematics of Computation,
%    Volume 64, Number 211, July 1995, pages 1147-1170.
%
%  Parameters:
%
%    Input, integer D, the spatial dimension.
%
%    Input, integer N, the maximum total degree.
%
%    Input, integer R, the number of monomials in D dimensions
%    of total degree N or less.
%
%    Input, integer ND, the number of data points.
%    It must be the case that ND <= R.
%
%    Input, real XD(D,ND), the data points, which must be distinct.
%
%    Output, integer PO(ND), the order (number of nonzero coefficients) for the
%    Lagrange basis polynomials.
%
%    Output, real PC(ND,R), the coefficients for the
%    Lagrange basis polynomials.
%
%    Output, integer PE(ND,R), the exponent indices for the
%    Lagrange basis polynomials.
%

%
%  Verify that R is correct.
%
  if ( r ~= mono_upto_enum ( d, n ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LAGRANGE_PARTIAL2 - Fatal error!\n' );
    fprintf ( 1, '  The value R is not correct.\n' );
    error ( 'LAGRANGE_PARTIAL2 - Fatal error!' );
  end

  if ( r < nd )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LAGRANGE_PARTIAL2 - Fatal error!\n' );
    fprintf ( 1, '  The value R = %d is less than ND = %d.\n', r, nd );
    error ( 'LAGRANGE_PARTIAL2 - Fatal error!' );
  end
%
%  Verify that the points are sufficiently distinct.
%
  [ d_min, d_max ] = r8col_separation ( d, nd, xd );
  d_tol = sqrt ( eps );

  if ( d_min < d_tol )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LAGRANGE_PARTIAL2 - Fatal error!\n' );
    fprintf ( 1, '  Some points are too close!\n' );
    fprintf ( 1, '  Minimum data point separation is = %g.\n', d_min );
    error ( 'LAGRANGE_PARTIAL2 - Fatal error!' );
  end
%
%  Initialize the polynomials Q to be all monomials of degree exactly N.
%
  qo = zeros ( r, 1 );
  qc = zeros ( r, r );
  qe = zeros ( r, r );

  for k = 1 : r
    qo(k,1) = 1;
    qc(k,1) = 1.0;
    qe(k,1) = k;
  end
%
%  Now set up the P polynomials.
%
  po = zeros ( nd, 1 );
  pc = zeros ( nd, r );
  pe = zeros ( nd, r );

  for k = 1 : nd
%
%  Find the polynomial Q(K:R)(X) which is most nonzero at X(K).
%
    i = r + 1;
    value_max = 0.0;

    for j = k : r

      o = qo(j);
      value = polynomial_value ( d, o, qc(j,1:o), qe(j,1:o), 1, xd(1:d,k) );

      if ( abs ( value_max ) <= abs ( value ) )
        i = j;
        value_max = value;
      end

    end

    if ( i == r + 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'LAGRANGE_PARTIAL2 - Fatal error!\n' );
      fprintf ( 1, '  I = R+1.\n' );
      fprintf ( 1, '  No polynomial was nonzero at the current node.\n' );
      error ( 'LAGRANGE_PARTIAL2 - Fatal error!' );
    end

    value = value_max;
%
%  Define P(K)(X) = Q(I)(X) / Q(I)(X(k)
%
    o = qo(i);
    po(k) = qo(i);
    pc(k,1:o) = qc(i,1:o) / value;
    pe(k,1:o) = qe(i,1:o);
%
%  Modify P(1:k-1)(X).
%
    for j = 1 : k - 1

      oj = po(j);
      ok = po(k);

      value = polynomial_value ( d, oj, pc(j,1:oj), pe(j,1:oj), 1, xd(1:d,k) );

      [ o, c, e ] = polynomial_axpy ( - value, ok, pc(k,1:ok), pe(k,1:ok), ...
        oj, pc(j,1:oj), pe(j,1:oj) );

      po(j) = o;
      pc(j,1:o) = c(1:o)';
      pe(j,1:o) = e(1:o)';

    end
%
%  Modify Q(I:downto:K+1)
%
    for j = i : -1 : k + 1

      oj = qo(j-1);
      ok = po(k);

      value = polynomial_value ( d, oj, qc(j-1,1:oj), qe(j-1,1:oj), ...
        1, xd(1:d,k) );

      [ o, c, e ] = polynomial_axpy ( - value, ok, pc(k,1:ok), pe(k,1:ok), ...
        oj, qc(j-1,1:oj), qe(j-1,1:oj) );

      qo(j) = o;
      qc(j,1:o) = c(1:o)';
      qe(j,1:o) = e(1:o)';

    end
%
%  Modify Q(I+1:R)
%
    for j = i + 1 : r

      oj = qo(j);
      ok = po(k);

      value = polynomial_value ( d, oj, qc(j,1:oj), qe(j,1:oj), ...
        1, xd(1:d,k) );

      [ o, c, e ] = polynomial_axpy ( - value, ok, pc(k,1:ok), pe(k,1:ok), ...
        oj, qc(j,1:oj), qe(j,1:oj) );

      qo(j) = o;
      qc(j,1:o) = c(1:o)';
      qe(j,1:o) = e(1:o)';

    end

  end
%
%  Get rid of tiny coefficients.
%
  for i = 1 : nd
    oi = po(i);
    [ o, c, e ] = polynomial_compress ( ...
      po(i), pc(i,1:oi), pe(i,1:oi) );
    po(i) = o;
    pc(i,1:o) = c(1:o)';
    pe(i,1:o) = e(1:o)';
  end

  return
end


text:
classdef (Abstract) Directive < handle

    properties(SetAccess = protected)
        Values
    end

    properties(Abstract, SetAccess = immutable)
        Name
    end

    methods(Abstract, Access = protected)
        [success, toBeGarbadeCollected] = tryToupdateFromValue(obj, value)
    end

    methods(Sealed)
        function [success, toBeGarbadeCollected] = updateFromValue(obj, value)
            [success, toBeGarbadeCollected] = obj.tryToupdateFromValue(value);
            if ~success
                warning('MBeautifier:Directive:InvalidValue', ['The value "', char(value), '" is not a valid directive value for directive "', obj.Name, '".']);
            end
        end
    end
end

author: davidvarga

repo: MBeautifier


text:
function [nmps2] = mGal2nmps2(mGal)
% Convert acceleration from milligalileos to nanometers per second squared
% Chad A. Greene 2012
nmps2 = mGal*1e+4;
end

author: Sable

repo: mcbench-benchmarks


text:
function y = cvx_constant( x )
narginchk(1,1);
y = x;

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [out_vec] = ch_model_tgn(in_vec)
% ch_model_tgn.m - TGn channel model simulation (load external coeff's)
global numPkts  PER_snr  snr_idx  ch_type  num_taps ...
       H  H_idx_inc  H_idx_ini  H_idx  H_cnt  H_dline;

% PER graph settings (actual values are user-selected (GUI))
%%%%% numPkts = 200;
%%%%% PER_snr = 29:2:37;

% To select channel matrices
numCoherenceTimes = ceil(10*numPkts/442);      % num of coherence times to test
numSampPerCohere = numPkts/numCoherenceTimes;  % max value = 442 (num of ch. samples in one coherence time)

% Input parameters
Clk = in_vec(1);
Ntx = in_vec(2);
Nrx = in_vec(3);
Nltf = in_vec(4);

% Index for additive noise
%     - add noise after 8 legacy sym, 2 HT-LTF symbols (with 6 sample delay)
%     - done to simulate with perfect channel estimation...
nse_idx = 1; %%5*80*(8+Nltf) + 5*6;

% Channel inputs
in_len = length(in_vec(5:end))/4;
ch_in  = zeros(4, in_len);
ch_in(1,:) = in_vec(0*in_len+3+(1:in_len)).';
ch_in(2,:) = in_vec(1*in_len+3+(1:in_len)).';
ch_in(3,:) = in_vec(2*in_len+3+(1:in_len)).';
ch_in(4,:) = in_vec(3*in_len+3+(1:in_len)).';
ch_in = ch_in(1:Ntx,:);  % keep only used Tx ant. inputs

% Channel outputs
ch_out = zeros(4, in_len);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Channel Model (external coeff's)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Various channel models, widths
ChD_Width   = [   10     10     10     10     10     10     10     10     10     20     30     30     30     40     50     50     50     50];

% Select channel model
if (1)
    ch_width = ChD_Width;
end

% Initialize channel model
if (Clk==0)
    % Initialize PER test
    snr_idx = 1;

    % Load Channel coeff's
    load ..\..\tgn_ch_taps\tgn_chD_nlos(fluor)_2x2_short;
    %load ..\..\tgn_ch_taps\tgn_chD_nlos_15m_4x4;

    % Setup channel coeff's for PER test
    H_cohere_prd = ceil(size(H, 4)/100);  % ch. mat samples in one coherence time (100 simulated)
    H_idx_inc = round(H_cohere_prd / numSampPerCohere);

    % Setup indices for channel coeff's
    num_taps = size(H, 3);
    H_idx_ini = 100;   % initial index for channel matrices
    H_idx = H_idx_ini;
    H_cnt = 1;

    % Initialize channel delay line
    tot_dly = sum(ch_width)/10;
    H_dline = zeros(Ntx, tot_dly);
end

% Use fixed channel for each packet
%     - assume 38400 samples for each packet (with silence added)
if 1  %(FIXED_CHANNEL)
 % Get channel coeff's
 H_tap(1:Ntx,1:Nrx,:) = H(1:Ntx,1:Nrx,:,H_idx);
 %H_tap(:,:,:) = H(:,:,:,H_idx);
 % Increment index after every packet...
 H_idx = H_idx + H_idx_inc;
 if (H_idx > size(H, 4)) H_idx = H_idx_ini; end

 % Compute channel power...
 H_pow = 0;
 for m=1:num_taps
    for k = 1:(ch_width(m)/10)
        He = H_tap(:,:,m);
        H_pow = H_pow + sum(sum(He.*conj(He)))/Nrx;
    end
 end
end

% Process all input samples...
for i=1:in_len

    % Add new sample to delay line
    H_dline = [ch_in(:,i) H_dline(:,1:end-1)];

    % Apply channel filter
    H_dline_idx = 1;
    H_out = zeros(Nrx,1);
    for m=1:num_taps
        % Apply tap to inputs (considering tap width)
        H_samp = zeros(Ntx, 1);
        for k = 1:(ch_width(m)/10)
            % 10 ns (100 MHz) sampling
            H_samp = H_samp + H_dline(:,H_dline_idx);
            H_dline_idx = H_dline_idx + 1;
        end
        H_out = H_out + H_tap(:,:,m).'*H_samp;
    end

    % Store channel output
    ch_out(1:Nrx,i) = H_out(1:Nrx);
end

% Compute noise variance (time-domain)
nse_sigma = H_pow*(1/5);                % signal upsampled by 5
nse_sigma = sqrt(0.5*nse_sigma);

% Add AWGN to Rx streams
snr_val = PER_snr(snr_idx);
for i=1:Nrx
     % Generate additive noise
     nse_size = size(ch_out(i,nse_idx:end));
     nse_out = nse_sigma * (randn(nse_size) + j*randn(nse_size));
     nse_out = nse_out * 10^(-snr_val/20);
     % Add noise to channel output
     ch_out(i,nse_idx:end) = ch_out(i,nse_idx:end) + nse_out;
end

% Next SNR Value for PER Curve...
if (mod(Clk,numPkts)==(numPkts-1) && snr_idx<length(PER_snr))
    snr_idx = snr_idx+1;
end

% Return output
out_vec = reshape(ch_out.', 4*in_len, 1);


text:
function  [fy,f]=FFT(y,Fs)

%==========================================================================
%  (Updated Version 03/09/2013)
%  Usage :
%  Function  [fy]=FFT(y,Fs)
%
%  1)computes the Power spectral density and Amplitude spectrum (P(f),F(f))
%  of 1d signal y(t) with sample rate  Fs (Nyquist rate) which is known
%  apriori. The results are plotted in 3 figures which correspond to simple
%  PSD,logarithmic PSD (dB) and  Amplitude Specturm respectively.
%                             _______________
%  such that  Ampitude(f) = \/    PSD(f)
%
%  2)The usefulness of this function is the adjustment of the frequency axis
%
%  3)The fast Fourier transform is computed with Matlab built-in function
%  fft, but for signals whose lengths <1000 points, one can use the nested
%  function   y=Fast_Fourier_Transform(X,N) .
%
%  Demo :
%
%                 Fs=800;
%                 Tf=2;
%                 t=0:1/Fs:Tf;
%                 f=[40 75];
%                 Amp=[4.5 9.22];
%                 sigma=1.33;
%                 y=Amp(1)*exp(j*2*pi*t*f(1))+Amp(2)*exp(j*2*pi*t*f(2));
%                 N=(sigma/sqrt(2))*(randn(size(t))+j*randn(size(t)));
%                 y=y+N;
%                 figure, plot(t,y),xlabel('time (s)'),ylabel('Voltage (v)'),
%                 title(strcat('Signal corrupted with AWGN, \sigma=',num2str(sigma))),
%                 fy=FFT(y,Fs);
%
% (c) KHMOU Youssef , Signal Processing 2013
%==========================================================================

if nargin<2
    error(' Sampling frequency is required to sompute (PSD,F(y)) ! ');
end

% In case that the input vector is matrix :  Maping with vect{} .
y=y(:).';
L=length(y);

%  (2^N) :Number of points for computing the FFT
N=ceil(log2(length(y)));

% FFT
fy=fft(y,2^N)/(L/2);
%------------------------------------------------
% for length<1000 one can replace fft with function :
% fy=Fast_Fourier_Transform(y,2^N)/(L/2); (line 84)
%------------------------------------------------

% Amplitude adjustment by checking for complex input y
if isreal(y)==0
    fy=fy/2;
end

% PSD
Power=fy.*conj(fy);
%Phase Angle
phy=angle(fy);

%  Frequency axis
f=(Fs/2^N)*(0:2^(N-1)-1);

%if nargin==4

% Figures------------------------------------------------------------------
ff1=figure;
plot(f,Power(1:2^(N-1)),'r'),  xlabel('  Frequency (Hz)'), ylabel(' Magnitude (w)'),
title('  Power Spectral Density'), grid on;
set(ff1,'Name','PSD');

ff2=figure;
plot(f,10*log10(Power(1:2^(N-1))),'r'),  xlabel('  Frequency (Hz)'), ylabel(' Magnitude  (dB)'),
title('  Power Spectral Density, logarithmic scale '), grid on;
set(ff2,'Name','10*log10(PSD)');

ff3=figure;
plot(f,sqrt(Power(1:2^(N-1))),'r'),  xlabel('  Frequency (Hz)'), ylabel('|F(Y)|'),
title('  Amplitude Spectrum'), grid on;
set(ff3,'Name','|F(y)|');

ff4=figure;
plot(f,phy(1:2^(N-1)),'b'), xlabel(' Frequency (Hz)'), ylabel(' arg(F(Y))'),
title(' Phase spectrum'), grid on;
set(ff4,'Name','arg(F(Y))');
%end
%==========================================================================
% function z=Fast_Fourier_Transform(x,nfft)
%
% N=length(x);
% z=zeros(1,nfft);
% Sum=0;
% for k=1:nfft
%     for jj=1:N
%         Sum=Sum+x(jj)*exp(-2*pi*j*(jj-1)*(k-1)/nfft);
%     end
% z(k)=Sum;
% Sum=0;% Reset
% end
% return
%=========================================================================

% Additional plot :  PDF
%M=max(real(y));
%M=2*M;
%x_range=-M:M/20:M;
%ff4=figure;
%hist(real(y),x_range),xlabel('Signal values'), ylabel('PDF');
%title('  Probability density function of the input signal'), grid minor;
%set(ff4,'Name','PDF');
%clear M x_range


text:
function varargout = plotlattice(varargin)
%PLOTLATTICE  Plots an integer lattice
%
% p = plotlattice(C,which,c,size,options)
%
% C    :  Constraint object
% which:  'inner' (default) or 'outer'
% color:  color [double] ([r g b] format) or char from 'rymcgbk'
% size :  Size of marker
% options: options structure from sdpsettings
% Example
% sdpvar x1 x2
% plot(x1^2+x2^2 <= 1.5,'green');hold on
% plotlattice(x1^2+x2^2 <= 1.5,'outer','yellow');
% plotlattice(x1^2+x2^2 <= 1.5,'inner','black');

h = ishold;
hold on
F = varargin{1};
if nargin > 1 && ~isempty(varargin{2})
    which = varargin{2};
else
    which = 'inner';
end
if nargin > 2 && ~isempty(varargin{3})
    color = varargin{3};
else
    color = 'yellow';
end
if nargin > 3 && ~isempty(varargin{4})
    size = varargin{4};
else
    size = 5;
end
if nargin > 4
    ops = varargin{5};
    if ~isempty(ops)
        ops = sdpsettings(ops,'verbose',0);
    else
        ops = sdpsettings('relax',2,'verbose',0);
    end
else
    ops = sdpsettings('relax',2,'verbose',0);
end

X = recover(depends(F));
[~,L,U] = boundingbox(F,ops);

x = floor(L(1)):ceil(U(1));
y = floor(L(2)):ceil(U(2));
ThreeD = length(L)>2;
if ThreeD
    z = floor(L(3)):ceil(U(3));
else
    z = 1;
end
V = [];
for i = x
    for j = y
        for k = z
        switch which
            case 'outer'
                if ThreeD
                    l = plot3(i,j,k,'or','MarkerSize',size);
                else
                    l = plot(i,j,'or','MarkerSize',size);
                end
                set(l,'MarkerFaceColor',color);
                V = [V [i;j;k]];
            case 'inner'
                if ThreeD
                    assign(X,[i;j;k]);
                else
                    assign(X,[i;j]);
                end
                p = check(F);
                if min(p) >= 0
                    if ThreeD
                        l = plot3(i,j,k,'or','MarkerSize',size,'MarkerFaceColor','yellow');
                        V = [V [i;j;k]];
                    else
                        l = plot(i,j,'or','MarkerSize',size,'MarkerFaceColor','yellow');
                        V = [V [i;j]];
                    end
                    set(l,'MarkerFaceColor',color);
                end
            otherwise
                error
        end
        end
    end
end
if ~h
    hold off;
end
if nargout > 0
    varargout{1} = V;
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function res = isstruct(this);

    res = 1;

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
contactList{1,1} = 'title 1'; contactList{1,2} = 'author name 1'; contactList{1,3} = 'spam@email.com';
contactList{2,1} = 'title 2'; contactList{2,2} = 'author name 2'; contactList{2,3} = 'spam@email.com';

subjectLine = 'title of the emails';
bodyLine = ['Dear %s,\n',...
'\n',...
'We will be organizing XXX\n',...
'Sincerely,\n',...
'The Organizers'];
cmdLine = 'mail -s "%s" %s < email.txt';

for itemID =1:size(contactList,1)
    paperTitle = contactList{itemID,1};
    authorName = contactList{itemID,2};
    emailAddress = contactList{itemID,3};

    % write the file
    fp = fopen('email.txt','w');
    fprintf(fp,bodyLine,authorName, paperTitle);
    fclose(fp);

    % run command
    cmd = sprintf(cmdLine,subjectLine,emailAddress);
    system(cmd);

    % delete file
    delete('email.txt');
end

author: jianxiongxiao

repo: ProfXkit


text:
function factors = ComputeAllSimilarityFactors (images, K)
% This function computes all of the similarity factors for the images in
% one word.
%
% Input:
%   images: An array of structs containing the 'img' value for each
%     character in the word.
%   K: The alphabet size (accessible in imageModel.K for the provided
%     imageModel).
%
% Output:
%   factors: Every similarity factor in the word. You should use
%     ComputeSimilarityFactor to compute these.
%
% Copyright (C) Daphne Koller, Stanford University, 2012

n = length(images);
nFactors = nchoosek (n, 2);

factors = repmat(struct('var', [], 'card', [], 'val', []), nFactors, 1);
k = 1;
% Your code here:
for i = 1:n-1
	for j = i+1:n
		factors(k) = ComputeSimilarityFactor(images,K,i,j);
		k+=1;
	end
end
end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function digits = isolateDigits(x, base)
digits = [];

if ~exist('base','var') || (base < 2 || base > 16)
    base = 10;
end
if base ~= uint32(base)
    return;
end
if x == 0
    digits = 0;
    return;
end

kk = 1;
while 1

    d = mod(x,base);
    x = floor(x/base);

    if (x == 0) && (d == 0)
        break;
    end

    digits(kk) = d;
    kk = kk+1;

end

author: BUNPC


text:
function y=stdtpdf(x,mu,sigma2,nu)
% Probability Density Function (PDF) for the Standardized T distribution
%
% USAGE:
%   Y = stdtpdf(X,MU,SIGMA2,NU)
%
% INPUTS:
%   X      - Standardized T random variables
%   MU     - Mean of X, either scalar or size(x)
%   SIGMA2 - Variance of X, either scalar or size(x)
%   NU     - Degree of freedom parameters, either scalar or size(x)
%
% OUTPUTS:
%   Y     - Probability density evaluated at X
%
% COMMENTS:
%   NU>2
%
% REFERENCES:
%   [1] Cassella and Berger (1990) 'Statistical Inference'
%
% See also STDTCDF, STDTINV, STDTRND, STDTLOGLIK, TPDF

% Copyright:
% Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 6    Date: 8/21/2014

[T,K]=size(x);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if K~=1
    error('X must be a column vector');
end

if nargin==4
    if length(mu)~=1 && ~all(size(mu)==[T K])
        error('mu must be either a scalar or the same size as X');
    end
    if any(sigma2<=0)
        error('sigma2 must contain only positive elements')
    end
    if length(sigma2)==1
        sigma2=sigma2*ones(T,K);
    elseif size(sigma2,1)~=T || size(sigma2,2)~=1
        error('sigma2 must be a scalar or a vector with the same dimensions as X');
    end
    if length(nu)>1 || nu<=2
        error('nu must be a scalar greater than 2');
    end
    x=x-mu;
else
    error('Only 4 inputs supported');
end

constant = exp(gammaln( 0.5 * (nu + 1)) - gammaln(0.5 * nu));
y = constant ./ sqrt(pi * (nu - 2) * sigma2) .* (1 + (x-mu) .^ 2.0 / (sigma2 * (nu - 2))) .^ (-(nu + 1) / 2);

author: bashtage

repo: mfe-toolbox


text:
%%  StoichTools: Tools for Doing Stoichiometry
%
% StoichTools comprises a set of Matlab functions for doing stoichiometric
% analysis. These functions parse standard chemical notation for a variety
% of stoichiometric calculations including finding molecular weights,
% balancing reactions for atom and charge conservation, finding independent
% reactions, and displaying formulas in Hill notation. The functions
% account for both change and atomic balances so they can be used to
% balance ionic reactions and chemical half reactions.
%
% StoichTools has extensive documentation including a set of worked
% homework problems demonstrating use of the functions.
%
% These functions were developed to support introductory courses in
% Chemical Engineering.
%
%  Jeff Kantor
%  December 18, 2010

%% What is StoichTools?
%
% StoichTools works with two types of data:
%
% # *Chemical formulas*. Each chemical formula is a string written in a
% nearly universal chemical notation. For example, |H2SO4| represents
% Sulfuric Acid. Grouping is allowed (e.g., |CH3(CH2)6CH3| for octane) with
% either parentheses '()' or brackets '[]'.  Charge is indicated by a
% trailing + or - followed by an optional number (e.g., |Fe+3| or |HSO4-|).
% Phase information may be included as a terminal (aq), (l), (g), or (s).
% Cell arrays can be used in most places to work with multiple formulas at
% one time (e.g., {'H2SO4','H+','SO4-2'}).
% # *Atomic representation*. Many calculations require knowledge of the
% charge, and of number of atoms of each type in a chemical species. This
% is maintained in a Matlab structure where r.C, for example, is the
% number of carbon atoms. The symbol after the dot is the standard 1 or 2
% character symbol for an element. The symbol |Q| is reserved to indicated
% charge. A Matlab structure array is used to store multiple atomic
% reprentations in a single variable.
%
% StoichTools provides functions for the following types of chemical
% calculations:
%
% *Working with Chemical Formulas*
%
% * |r = parse_formula(s)| processes a chemical formula to produce an
%   atomic representation. This function is mainly used by other functions
%   to process chemical formulas.
% * |hillformula| processes a chemical formula or atomic reprentation to
%   produce a chemical formula in standard Hill notation. The Hill notation
%   widely used to represent species in chemical databases, such as the
%   NIST Chemistry Webbook.
%
% *Calculating Molecular Weights*
%
% * |mw = molweight(s)| computes the molecular weights of chemical
%   compounds. Input can be a chemical formula, a cell array of chemical
%   formulas, or an array of atomic representations. If no output is
%   indicated, then a table of molecular weights is printed.
%
% *Stoichiometry*
%
% * |[A,atoms,species] = atomic(s)| constructs the atomic matrix for a set
%   of chemical compounds. Element |A(i,j)| is the number of |atoms{i}| in
%   |species{j}|. Inputs may be chemical formula, a cell array of chemical
%   formulas, If there are ionic species, then a special atom 'Q' is
%   indicates the charge of the species. If no output is indicated, then
%   the atomic matrix is displayed in tabular form.
% * |V = stoich(s)| computes the stoichiometric matrix for a set of
%   chemical compounds. The input is a cell array of chemical formulas, or
%   an array of atomic representations. The columns of |V| correspond to
%   independent chemical reactions satisfying atomic and charge balances.
%   Element |V(j,k)| is the stoichiometric coefficient for species |j| in
%   reaction |k|. A negative value denotes a reactant, a positive value
%   denotes a product. If no output is indicated, then |disp_reaction| is
%   used to display all independent reactions.
% * |Vout = disp_reaction(V,s)| If no output is indicated, then format
%   and displays the chemical reactions denoted by stoichiometric matrix
%   |V| and the array of species |s|. The species may be cell array of
%   formulas or an array of atomic representations. If feasible, the
%   coefficients are scaled to integers.  It integer coefficients are too
%   long, then either rational or floating point coefficients are
%   displayed. If an output is indicated, then |Vout| is a stoichiometric
%   matrix with rescaled coefficients, and the reactions are not displayed.
%
% *Homework Problems with Solutions*
%
% The StoichTools folder includes a number of worked homework problems.
% These are Matlab scripts with titles in the pattern |HW_xx.m|. Each
% script begins with a cell containing the problem statement. Subsequent
% cells demonstrate solution to the problem. The homework files can be
% sviewed by using the Matlab publishing function.

%% Parsing Chemical Formulas
%
% Given a set of chemical species, |r = parse_formula(s)| parses a cell
% array of chemical formulas to produce a structure array r. The value is
% the number of atoms of that element present in the corresponding formula.
% The structure array includes a field for each atomic element in the set
% of species. We call this the atomic represenation of the species.

% Parsing methane

parse_formula('CH4')

%% Additional Parsing Examples

ex{1} = 'NaHCO3';
ex{2} = 'KFe3(SO4)2(OH)6';     % Jorosite
ex{3} = 'KFe3(AsO4)2(HAsO4)2'; % Potassium-Iron-Arsenate
ex{4} = '(CH4)8(H2O)46';       % Methane Clathrate
ex{4} = 'HSO4-(aq)';

for k = 1:length(ex)
    disp(ex{k});
    parse_formula(ex{k})
end

%% Chemical Abbreviations and Isotopes
%
% * Formulas may include D (Deuterium) or T (Tritium). These are treated as
%   elements and included as distinct species in any atom balances.
% * The common organic chemistry abbreviations Me (Methyl, CH3), Et (Ethyl,
%   C2H5), Bu (Butyl, C4H9), Ph (Phenol, C6H5) may be included in formulas.
%   These are replaced by their atomic formulas during the parsing process.
% * The symbols M (any metal) and X (any halogen) may be used in formulas.
%   Formulas containing the symbol M or X have unknown molecular weight.

parse_formula('D2O')
parse_formula('EtOH')
molweight({'H2O','D2O','T2O','EtOH','PhOH','TiO2','MO2'});

%% Non-stoichiometric Formulas
%
% Some applications of stoichiometry involve complex chemical compounds not
% easily described by simple chemical fomulas. So-called
% 'non-stoichiometric' compounds can be also be parsed.

bacteria = 'CH1.8N0.24O0.36';
parse_formula(bacteria);

%% From Atoms to Chemical Formulas
%
% Given a structure array of atomic representations, |s = hillformula(r)}
% constructs a cell array of corresponding chemical formulas.

% Formula for octane

octane.C = 8;
octane.H = 18;
hillformula(octane)

%% Hill Notation & Canonical Representations
%
% The Hill notation is a commonly used system for writing chemical formulas
% in a standard form. % |hillformula(r)| produces a simple canonical
% representation of a chemical species. Note, however, that there may be
% many isomers for a given formula.

s = {'Zr3B2','HBr','HCl','CH3(CH2)6CH3','NaCO3','CaC2','CH3OH', ...
     'CH3COOH','HNO3','H2SO4','NH3','SnH4','CH3HgCH3','(CH3CH2)4Pb', ...
     '[Co(NH3)6]+3','[B12H12]-2'};

fprintf('\n%-15s %-15s\n----------      ----------\n', ...
    'Formula','Hill Notation');
for k = 1:length(s)
    fprintf('%-15s %-15s\n',s{k},char(hillformula(s{k})));
end

%% Molecular Weight
%
%  mw = molweight(s)
%  mw = molweight(r)
%
% Given a cell array of chemical formulas, or a structure array of atomic
% representations, |molweight| computes a corresponding vector of molecular
% weights.

% Molecular Mass of Dimethyl Mercury

s = 'CH3HgCH3';
mw = molweight('CH3HgCH3');
fprintf('Molecular Weight of Dimethyl Mercury (%s) = %g\n',s,mw);

%% Creating Molecular Weight Tables
%
% If molweight as no output, then it prints a table of molecular weights.

molweight(s);

%% Atomic Matrix
%
%  [A,atoms,species] = atomic(s)
%  [A,atoms,species] = atomic(r)
%
% Given a cell array of chemical formulas |s|, or a structure array of
% atomic representations |r|, |atomic| computes the atomic matrix A.
% |atoms| is a a cell array of the atomic elements, |species| is a cell
% array of species. A(i,j) is the number of atoms of element atoms{i} in
% species species{j}. % When called without an output argument, |atomic|
% displays the atomic matrix.

s = {'CH4','O2','H2O','CO2'};

atomic(s);
A = atomic(s);
disp(' ');
disp('A = ');
disp(A);

%% Atomic Matrix for Ionic Species
%
% For ionic species an additional row is added, labeled by 'Q', indicating
% the net charge on each of the species included in the matrix.

s = {'Fe+3','SO4-2','H+','OH-','H2O','Fe2(SO4)3'};
atomic(s);

%% Balancing a Reaction
%
% Given a cell array of chemical formulas, or an array of atomic
% representations, |stoich(s)| computes stoichiometric coefficients that
% satisfy charge and atom balances. If no output is specified, then
% balanced reactions are displayed.

stoich({'NaPb','CH3CH2Cl','(CH3CH2)4Pb','NaCl','Pb'});
stoich({'H+(aq)','OH-(aq)','H2O(l)'});

%% Stoichiometric Matrix
%
% Given a cell array of chemical formulas, or a structure array of atomic
% representations, |V = stoich(s)| computes the stoichiometric matrix |V|.
% |V(n,r)| is the stoichiometric coeffient of species |n| in reaction |r|.
% The atomic and stoichiometric matrices satisfies the relationship |A*V =
% 0|.

s = {'C8H18','O2','C','CO','CO2','H2O'};
V = stoich(s);
disp('Stoichiometric Matrix V = ');
disp(V);

%% Mulitple Independent Reactions
%
%  V = stoich(s)
%  disp_reaction(V,s)
%
% The columns of the stoichiometric matrix |V| represent independent
% reactions. The function |disp_reaction(V,s)| displays the reactions in a
% conventional human readable form.

s = {'C8H18','O2','C','CO','CO2','H2O'};
V = stoich(s);
disp_reaction(V,s);

%% Further Examples of Complex Reactions
%
% Examples from
% <http://www.chemistryhelp.net/chemistry-calculator/chemical-equation-balancer>

stoich({'P2I4','P4','H2O','H3PO4','PH4I'});

stoich({'[Cr(N2H4CO)6]4[Cr(CN)6]3','KMnO4','H2SO4','K2Cr2O7', ...
     'MnSO4','CO2','KNO3','K2SO4','H2O'});

stoich({'Cu(s)','HNO3(aq)','Cu(NO3)2(aq)','NO(g)','H2O(l)'});

stoich({'Cu','HNO3','H2O','Cu(NO3)2','NO'});

stoich({'KMnO4','C3H5(OH)3','K2CO3','Mn2O3','CO2','H2O'});

stoich({'K2Cr2O7','FeCl2','HCl','KCl', ...
    'CrCl3','FeCl3','H2O'});

stoich({'Bi(NO3)3(H2O)5','NaOH','H2O2','RuCl3', ...
     'NaNO3','NaCl','Bi2Ru2O7','H2O'});

stoich({'(NH4)2MoO4','NH4NO3','Na3PO4','H2O', ...
    '(NH4)3[P(Mo3O10)4]','NaNO3','NH3'});

stoich({'H2','Ca(CN)2','NaAlF4','FeSO4','MgSiO3','KI','H3PO4', ...
    'PbCrO4','BrCl','CF2Cl2','SO2','PbBr2','CrCl3','MgCO3', ...
    'KAl(OH)4','Fe(SCN)3','PI3','Na2SiO3','CaF2','H2O'});

stoich({'NH4ClO4','NaY(OH)4','Ru(SCN)3','PBr5','TiCl2CrI4','BeCO3', ...
    'Rb2ZrO3','ZnAt2','CAt2I2','Rb0.998YAt4','RuS2','BeZrO3','Zn(CN)2', ...
    'NaHBr1.997','H3PO4','TiCrO4','ClI','H2SO4','H2O'});

%% Chemical Equations with Ionic Charges
%
% The charge on ionic species is indicated by + or - followed by an
% optional digit indicating the amount of charge.  If ionic species are
% present, then a charge balance is include in the computation of the
% stoichiometric coefficients.

stoich({'ClO2+(aq)','H3O+(aq)','Cl2(g)','H2O(l)','ClO3-(aq)','ClO2(aq)'});
stoich({'Bi+3(aq)','HSnO2-(aq)','OH-(aq)','Bi(s)','H2O','SnO3-2(aq)'});
stoich({'CH3CH2OH','Cr2O7-2','H+','CH3COOH','Cr+3','H2O'});
stoich({'I-','I2','Mn+2','MnO4-','H+','H2O'});
stoich({'Cl2','Cl-','Fe+2','Fe+3'});
stoich({'Mn+2','BiO-3','H+','MnO4-','Bi3+','H2O'});
stoich({'NpO2+2','NpO2(OH)H2C2O4+','NpO2+','CO2','H+','O2'});
stoich({'H3PO4','(NH4)6Mo7O24','H+','(NH4)3PO4(MoO3)12','NH4+','H2O'});

%% Chemical Half Equations
%
% Include the bare electron 'e-' to balance chemical half reactions. In
% acidic solutions, if one of the main reactants contains oxygen, add 'H+'
% and 'H2O'. In basic solutions, if one of the main reactants contains
% oxygen then add 'OH-' and 'H2O'.

stoich({'Al+3(aq)','Al(s)','e-'});
stoich({'Cl-(aq)','Cl2(g)','e-'});

% Acidic Solutions

stoich({'MnO4-(aq)','Mn+2(aq)','H2O(l)','H+(aq)','e-'});
stoich({'O2(g)','H2O(l)','H+(aq)','e-'});
stoich({'Ag2O3','Ag+','H2O','H+','e-'});
stoich({'S2O3-2(aq)','S(s)','H2O(l)','H+(aq)','e-'});
stoich({'HOOCCOOH(aq)','CO2(g)','H2O(l)','H+(aq)','e-'});

% Alkali Solutions

stoich({'MnO4-(aq)','Mn+2(aq)','H2O(l)','OH-(aq)','e-'});
stoich({'Cr(OH)6-2','CrO4-2','H2O','OH-','e-'});
stoich({'NH3OH(aq)','N2(g)','H2O(l)','OH-(aq)','e-'});
stoich({'Al(OH)4-(aq)','Al(s)','H2O(l)','OH-(aq)','e-'});
stoich({'ZrO(OH)2','Zr','H2O','OH-','e-'});

%% Nested Formulas
%
% Matlab regular expressions capabilities are used to parse chemical
% formulas. While this keeps StoichTools simple and fast, one of the
% drawbacks of regular expressions is the difficulty of matching nested
% expressions. Thus nesting is limited to bracketed expressions inside of
% parentheses, or parentheses inside of brackets. By this rule, [Fe2(SO4)3]
% and (Fe2[SO4]3) are allowed, but (Fe2(SO4)3) and [Fe2[SO4]3] are not. In
% practice, chemical formula rarely need more than two levels of nesting.

disp('These work fine.');
molweight({'[Fe2(SO4)3]','(Fe2[SO4]3)'});

fprintf('\n\n');
try
    molweight({'(Fe2(SO4)3)','[Fe2[SO4]3]'})
catch exception
    disp('But this does not.');
    disp(exception.message);
end

%% Version History
%
% * 2010/12/18  Submitted to Matlab Central
% * 2010/12/19  Updated documentation, added solved homeworks
% * 2010/12/19  Put rows of the atomic matrix in Hill order
% * 2010/12/19  Expanded regular expression parsing to include phases
% * 2010/12/20  Enhanced parser to accept non-stoichiometric formulas
% * 2010/12/20  Enhanced disp_reaction for better coefficient formatting
% * 2010/12/21  Parser to include common symbols D, T, Et, Me, Bu, Ph
% * 2010/12/30  Fixed all mlint messages, reduced McCabe complexity
% * 2010/12/30  Update to Matlab Central
% * 2010/12/30  Further improvements to error handling (assert's)
% * 2010/12/31  Fixed bug with NaN in molweight
% * 2010/12/31  Renamed homework files so it makes more sense on MC
% * 2010/12/31  Update to Matlab Central
%
%
% To Do's
%
% * Add Generation/Consumption Analysis
% * Add Extent of Reaction Analysis
% * Include an electrochemistry howework example (battery?)
% * Add a display feature for stoich
% * Add webbook lookup for chemical property data

author: opencobra

repo: cobratoolbox


text:
function [Fp, Fi] = tt(F)
  warning('tt is deprecated, use triangle_triangle_adjacency instead.');
  [Fp, Fi] = triangle_triangle_adjacency(F);
end

author: alecjacobson

repo: gptoolbox


text:
function  fa = cmd3(a, om, omt)
% calculate the growth factor using 1/(ah)^3
% t0 is (1 - om)/om with om the relative matter density

%fa = (a ./ (1 + a .^3 * t0)) .^ 1.5;

fa = (a ./ (om + (omt - om) .* a.^3  + (1 - omt) .* a))  .^1.5;


text:
% STD_CLUSTMAXELEC - function to find the electrode with maximum absolute projection
%                  for each component of a cluster
% Usage:
%     >> [eleclist, setlist, complist] = std_clustmaxelec(STUDY, ALLEEG, clustind);
%
% Inputs:
%   STUDY        - STUDY set structure containing (loaded) EEG dataset structures
%   ALLEEG       - ALLEEG vector of EEG structures, else a single EEG dataset.
%   clustind     - (single) cluster index
%
% Outputs:
%   eleclist     - [cell] electrode list
%   setlist      - [integer] set indices for the cluster
%   complist     - [integer] component indices for the cluster
%
% Authors: Claire Braboszcz & Arnaud Delorme , CERCO, UPS/CRNS, 2011

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

function std_clustmaxelec(STUDY, ALLEEG, clusterind);

if nargin < 1
    help findicmaxelec;
    return;
end

fprintf('Finding electrodes with max weight for cluster %d\n', clusterind);
fprintf('-------------------------------------------------\n');

for index = 1:length(STUDY.cluster(clusterind).comps)
    set  = STUDY.cluster(clusterind).sets(1,index);
    comp = STUDY.cluster(clusterind).comps( index);

    [tmp maxelec] = max( abs(ALLEEG(set).icawinv(:, comp)) );
    indelec = ALLEEG(set).icachansind(maxelec);
    maxallelec{index} = ALLEEG(set).chanlocs(indelec).labels;
    allelec = unique_bc(maxallelec);

    fprintf('The electrode with the max weight for component %d of dataset %d is "%s"\n', comp, set, maxallelec{index});
end

for indelec=1:length(allelec)
    nbelec{indelec} = length(find(strcmp(allelec{indelec}, maxallelec) == 1));

    fprintf('Number of occurrence of electrode %s: %d\n', allelec{indelec}, nbelec{indelec});
end


text:
%DEMOLGD1 Demonstrate simple MLP optimisation with on-line gradient descent
%
%	Description
%	The problem consists of one input variable X and one target variable
%	T with data generated by sampling X at equal intervals and then
%	generating target data by computing SIN(2*PI*X) and adding Gaussian
%	noise. A 2-layer network with linear outputs is trained by minimizing
%	a  sum-of-squares error function using on-line gradient descent.
%
%	See also
%	DEMMLP1, OLGD
%

%	Copyright (c) Ian T Nabney (1996-2001)

% Generate the matrix of inputs x and targets t.

ndata = 20;			% Number of data points.
noise = 0.2;			% Standard deviation of noise distribution.
x = [0:1/(ndata - 1):1]';
randn('state', 42);
rand('state', 42);
t = sin(2*pi*x) + noise*randn(ndata, 1);

clc
disp('This demonstration illustrates the use of the on-line gradient')
disp('descent algorithm to train a Multi-Layer Perceptron network for')
disp('regression problems.  It is intended to illustrate the drawbacks')
disp('of this algorithm compared to more powerful non-linear optimisation')
disp('algorithms, such as conjugate gradients.')
disp(' ')
disp('First we generate the data from a noisy sine function and construct')
disp('the network.')
disp(' ')
disp('Press any key to continue.')
pause
% Set up network parameters.
nin = 1;			% Number of inputs.
nhidden = 3;			% Number of hidden units.
nout = 1;			% Number of outputs.
alpha = 0.01;			% Coefficient of weight-decay prior.

% Create and initialize network weight vector.
net = mlp(nin, nhidden, nout, 'linear');
% Initialise weights reasonably close to 0
net = mlpinit(net, 10);

% Set up vector of options for the optimiser.
options = foptions;
options(1) = 1;			% This provides display of error values.
options(14) = 20;		% Number of training cycles.
options(18) = 0.1;		% Learning rate
%options(17) = 0.4;		% Momentum
options(17) = 0.4;		% Momentum
options(5) = 1; 		% Do randomise pattern order
clc
disp('Then we set the options for the training algorithm.')
disp(['In the first phase of training, which lasts for ',...
    num2str(options(14)), ' cycles,'])
disp(['the learning rate is ', num2str(options(18)), ...
    ' and the momentum is ', num2str(options(17)), '.'])
disp('The error values are displayed at the end of each pass through the')
disp('entire pattern set.')
disp(' ')
disp('Press any key to continue.')
pause

% Train using online gradient descent
[net, options] = olgd(net, options, x, t);

% Now allow learning rate to decay and remove momentum
options(2) = 0;
options(3) = 0;
options(17) = 0.4;	% Turn off momentum
options(5) = 1;		% Randomise pattern order
options(6) = 1;		% Set learning rate decay on
options(14) = 200;
options(18) = 0.1;	% Initial learning rate

disp(['In the second phase of training, which lasts for up to ',...
    num2str(options(14)), ' cycles,'])
disp(['the learning rate starts at ', num2str(options(18)), ...
    ', decaying at 1/t and the momentum is ', num2str(options(17)), '.'])
disp(' ')
disp('Press any key to continue.')
pause
[net, options] = olgd(net, options, x, t);

clc
disp('Now we plot the data, underlying function, and network outputs')
disp('on a single graph to compare the results.')
disp(' ')
disp('Press any key to continue.')
pause

% Plot the data, the original function, and the trained network function.
plotvals = [0:0.01:1]';
y = mlpfwd(net, plotvals);
fh1 = figure;
plot(x, t, 'ob')
hold on
axis([0 1 -1.5 1.5])
fplot('sin(2*pi*x)', [0 1], '--g')
plot(plotvals, y, '-r')
legend('data', 'function', 'network');
hold off

disp('Note the very poor fit to the data: this should be compared with')
disp('the results obtained in demmlp1.')
disp(' ')
disp('Press any key to exit.')
pause
close(fh1);
clear all;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function R = functionRlocalscattering3D(M_H, M_V, d_H, d_V, varphi, stdphi, theta, stdtheta, dist)
%Generate the spatial correlation matrix for the 3D local scattering model,
%defined in (7.17) for different angular distributions.
%
%INPUT:
% M_H           = Number of antennas per horizontal row, i.e., # of columns
% M_V           = Number of horizontal rows
% d_H           = Horizontal antenna spacing in multiples of wavelength
% d_V           = Vertical antenna spacing in multiples of wavelength
% varphi        = Azimuth angle in radians
% stdphi        = Azimuth angular spread (standard deviation or limit)
% theta         = Elevation angle in radians
% stdtheta      = Elevation angular spread (standard deviation or limit)
% dist          = 'Gaussian' or 'Laplace' or 'Uniform'
%
%OUTPUT:
% R             = M x M channel covariance matrix, where M = M_H x M_V
%
%
%This Matlab function was developed to generate simulation results to:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

% Index mapping functions
i = @(m) mod(m-1, M_H);
j = @(m) floor((m-1)/M_H);

%Create the correlation matrix
M = M_H*M_V;
R = zeros(M);

% Define angle distributions as functions
if (strcmp(dist,'Laplace'))
    f1 = @(x) exp(-sqrt(2)*abs(x-varphi)/stdphi)/(sqrt(2)*stdphi);
    f2 = @(x) exp(-sqrt(2)*abs(x-theta)/stdtheta)/(sqrt(2)*stdtheta);
elseif (strcmp(dist,'Gaussian'))
    f1 = @(x) exp(-(x-varphi).^2/(2*stdphi^2))/(sqrt(2*pi)*stdphi);
    f2 = @(x) exp(-(x-theta).^2/(2*stdtheta^2))/(sqrt(2*pi)*stdtheta);
elseif (strcmp(dist,'Uniform'))
    f1 =  @(x)1/(2*stdphi);
    f2 =  @(x) 1/(2*stdtheta);
else
    error('Please provide a valid angle distribution')
end

%% Compute correlation matrix
LOOKUP = zeros(2*M_H-1,M_V); %Define a lookup table to speedup the computation process
for m = 1:M
    for l = 1:m
        indi = i(l)-i(m) + M_H;
        indj = j(l)-j(m) + M_V;
        if (LOOKUP(indi,indj)==0)
            integrand = @(p,t) exp(1j*2*pi*d_V*(j(l)-j(m)).*sin(t)) .* exp(1j*2*pi*d_H*(i(l)-i(m)).*cos(t).*sin(p)) .* f1(p) .* f2(t);
            if (strcmp(dist,'Uniform'))
                LOOKUP(indi,indj) = integral2(integrand, varphi-stdphi, varphi+stdphi, theta-stdtheta, theta+stdtheta);
            else
                LOOKUP(indi,indj) = integral2(integrand, varphi-20*stdphi, varphi+20*stdphi, theta-20*stdtheta, theta+20*stdtheta);
            end
        end
        R(m,l) = LOOKUP(indi,indj);
        R(l,m) = R(m,l)';
    end
end

author: emilbjornson

repo: massivemimobook


text:
function out = spm_dartel_norm(job)
% Warp individuals to match template
% FORMAT spm_dartel_norm(job)
% job.flowfields - Flow-fields
% job.images     - Image to warp
% job.interp     - Interpolation method
% job.K          - 2^K timesteps are used
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_dartel_norm.m 5668 2013-10-03 18:34:18Z guillaume $

PU = job.flowfields;
PI = job.images;
jactransf = job.jactransf;
intrp    = job.interp;
K        = job.K;

for i=1:numel(PU),
    drawnow;
    NU = nifti(PU{i});
    [pth,nam,ext,num] = spm_fileparts(NU.dat.fname);
    fprintf('%s: ',nam);
    if jactransf,
        [y,dt] = spm_dartel_integrate(NU.dat,[1 0], K);
        dt     = max(dt,eps);
        if jactransf~=1,
            dt = dt.^jactransf;
        end
    else
        y      = spm_dartel_integrate(NU.dat,[1 0], K);
    end;
    y1 = double(y(:,:,:,1));
    y2 = double(y(:,:,:,2));
    y3 = double(y(:,:,:,3));

    for m=1:numel(PI),
        [pth1,nam,ext,num] = spm_fileparts(PI{m}{i});
        NI = nifti(fullfile(pth1,[nam ext]));
        NO = NI;
        if jactransf,
            NO.dat.fname=fullfile(pth,['mw' nam ext]);
            NO.dat.scl_slope = 1.0;
            NO.dat.scl_inter = 0.0;
            NO.dat.dtype     = 'float32-le';
        else
            NO.dat.fname=fullfile(pth,['w' nam ext]);
        end;
        NO.dat.dim = [NU.dat.dim(1:3) NI.dat.dim(4:end)];
        NO.mat  = NU.mat;
        NO.mat0 = NU.mat;
        NO.mat_intent  = 'Aligned';
        NO.mat0_intent = 'Aligned';
        NO.descrip = 'Dartel warped';
        NO.extras  = [];
        create(NO);
        fprintf('%s',nam); drawnow;

        for j=1:size(NI.dat,4),
            if sum(sum((NI.mat  - NU.mat ).^2)) < 0.0001 && ...
               sum(sum((NI.mat0 - NU.mat0).^2)) < 0.0001,
                ty1 = y1;
                ty2 = y2;
                ty3 = y3;
            else
                mat = NI.mat;
                if ~isempty(NI.extras) && isstruct(NI.extras) && isfield(NI.extras,'mat'),
                    mat1 = NI.extras.mat;
                    if size(mat1,3) >= j && sum(sum(mat1(:,:,j).^2)) ~=0,
                        mat = mat1(:,:,j);
                    end;
                end;
                M   = mat\NU.mat0;
                ty1 = M(1,1)*y1 + M(1,2)*y2 + M(1,3)*y3 + M(1,4);
                ty2 = M(2,1)*y1 + M(2,2)*y2 + M(2,3)*y3 + M(2,4);
                ty3 = M(3,1)*y1 + M(3,2)*y2 + M(3,3)*y3 + M(3,4);
            end;
            for k=1:size(NI.dat,5),
                for l=1:size(NI.dat,6),
                    f  = NI.dat(:,:,:,j,k,l);
                    spl_param = [intrp,intrp,intrp,0,0,0];
                    if intrp>1,
                        f = spm_bsplinc(f,spl_param);
                    end;
                    f = spm_bsplins(f,ty1,ty2,ty3,spl_param);
                    if jactransf,
                        NO.dat(:,:,:,j,k,l)=f.*dt;
                    else
                        NO.dat(:,:,:,j,k,l)=f;
                    end;
                    fprintf('\t%d,%d,%d', j,k,l); drawnow;
                end;
            end;
        end;
        fprintf('\n'); drawnow;
    end;
end;

PU = job.flowfields;
PI = job.images;
jactransf = job.jactransf;
out.files = cell(numel(PU),numel(PI));
for i=1:numel(PU),
    [pth,nam] = spm_fileparts(PU{i});
    for m=1:numel(PI),
        [pth1,nam,ext,num] = spm_fileparts(PI{m}{i});
        if jactransf,
            fname = fullfile(pth,['mw' nam ext]);
        else
            fname = fullfile(pth,['w' nam ext]);
        end;
        out.files{i,m} = fname;
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%Needed linguistic variables for importance weight of each criterion as generalized fuzzy
%numbers:
VL=[0 0 0.1 0.2 1];     %Very low
L=[0.1 0.2 0.2 0.3 1];  %Low
ML=[0.2 0.3 0.4 0.5 1]; %Medium low
M=[0.4 0.5 0.5 0.6 1];  %Medium
MH=[0.5 0.6 0.7 0.8 1]; %Medium high
H=[0.7 0.8 0.8 0.9 1];  %High
VH=[0.8 0.9 1.0 1.0 1]; %Very High

%Needed linguistic variables for ratings as fuzzy numbers:

VP=[0,0,1,2, 1];     % Very Poor
P=[1,2,2,3, 1];     % Poor
MP=[2,3,4,5, 1];     % Medium Poor
F=[4,5,5,6, 1];     % Fair
MG=[5,6,7,8, 1];    % Medium Good
G=[7,8,8,9, 1];     % Good
VG=[8,9,10,10,1];   % Very Good


text:
function [ x, odata, opts ] = solver_TraceLS( A, b, lambda, x0, opts )
% SOLVER_TRACELS Unconstrained form of trace-regularized least-squares problem.
% [ x, odata, opts ] = solver_TraceLS( A, b, lambda, x0, opts )
%    Solves the trace-regularized least squares problem
%        minimize (1/2)*norm( A * X - b )^2 + lambda * trace( X )
%        with the constraint that X is positive semi-definite.
%    A must be a matrix or a linear operator, b must be a vector,
%    and lambda must be a real positive scalar.
%    The initial point x0 and option structure opts are
%    both optional.
%
%    If "A" is a sparse matrix, and nnz(Z) = length(b), then it is assumed
%    that the nonzero entries correspond to samples, and then the
%    corresponding sampling operator is used.
%
%    If opts.largescale = true, then uses an iterative method
%    to compute the eigenvalue decomposition used with prox_trace.
%
%   See also solver_L1RLS (aka The Lasso)

% Added Feb 7, 2011
error(nargchk(3,5,nargin));
if nargin < 4, x0 = []; end
if nargin < 5, opts = []; end
if ~isfield( opts, 'restart' ),
    opts.restart = 100;
end
if issparse(A) && nnz(A) == length(b)
    A = linop_subsample(A);
end
if isfield( opts,'largescale')
    prx = prox_trace( lambda, opts.largescale );
    opts = rmfield(opts,'largescale');
else
    prx = prox_trace( lambda );
end
% Note: the proximity operator of trace is really a combination
%  of the proximity operator of trace and the indicator
%  set of positive semi-definite matrices.
%  So we do not need to *explicitly* include the positive semi-definite constraint.

[x,odata,opts] = tfocs( smooth_quad, { A, -b }, prx, x0, opts );

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function d = dot(m1,m2,varargin)
% inner product between two Miller indece
%
% Syntax
%   d = dot(m1,m2)
%   d = dot(m1,m2,'antipodal')
%
% Input
%  m1, m2 - @Miller
%
% Output
%  d - double, same size as |m1| and |m2|
%
% Options
%  noSymmetry - do *not* consider sym. equiv. directions
%  max       - (default) maximum dot product with respect to all sym. equiv.
%  min       - minimum dot product with respect to all sym. equiv.
%  all       - all dot products with respect to sym. equiv.
%  antipodal - include antipodal symmetry
%

% maybe we should ignore symmetry
if check_option(varargin,'noSymmetry') || ~isa(m2,'Miller')
  d = dot@vector3d(m1,m2,varargin{:});
  return
end

% if we should consider symmetry - it must be the same on both sides
if m1.CS ~= m2.CS, warning('Symmetry mismatch'); end

% maybe we should return a full matrix of dot products to all symmetrically
% equivalent directions
if check_option(varargin,'all')

  if length(m1) == 1
    m1 = repmat(m1,size(m2));
  else
    m2 = repmat(m2,size(m1));
  end

elseif (length(m1)==1 || length(m2) == 1) % use dot_outer whenever possible
  d = dot_outer(m1,m2,varargin{:});

  if length(m1) == 1
    d = reshape(d,[size(m2),size(d,3)]);
  else
    d = reshape(d,[size(m1),size(d,3)]);
  end
  return
end

% symmetrize
s = size(m1);
m1 = symmetrise(m1,varargin{:});
m2 = repmat(reshape(m2,1,[]),size(m1,1),1);

% vector3d dot product
d = dot@vector3d(m1,m2,varargin{:});

% which angle to return
if check_option(varargin,'min')
  d = reshape(min(d,[],1),s); % minimum angle of all symmetricaly equ.
elseif ~check_option(varargin,'all')
  d = reshape(max(d,[],1),s); % maximum angle of all symmetricaly equ.
end

author: mtex-toolbox

repo: mtex

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


text:
%maps label strings to integers
classdef Labeling < handle

    properties (Access = public, Constant)
        kIgnoreClass = -2;
        kInvalidClass = -1;
        kNullClass = 0;
        kNullClassStr = 'NULL';
        kIgnoreStr = 'ignore';
    end

    properties (GetAccess = public)
        numClasses;
        classNames;
    end

    properties (Access = private)
        classesMap;
    end

    methods (Access = public)
        function obj = Labeling(classNames)
            obj.classNames = classNames;
            if ~isempty(classNames)
                obj.numClasses = length(classNames);
                obj.createLabeling(classNames);
            end
        end

        function valid = isValidLabel(obj,labelStr)
            valid = isKey(obj.classesMap,labelStr);
        end

        function classStr = stringForClassAtIdx(obj,idx)
            if idx == Labeling.kNullClass
                classStr = Labeling.kNullClassStr;
            elseif idx == Labeling.kIgnoreClass
                classStr = Labeling.kIgnoreStr;
            else
                classStr = obj.classNames{idx};
            end
        end

        function idx = idxOfClassWithString(obj,classStr)
            if isempty(obj.classesMap)
                idx = [];
            else
                if ~isKey(obj.classesMap,classStr)
                    fprintf('%s: %s\n',Constants.kUndefinedClassError,classStr);
                    idx = [];
                else
                    idx = obj.classesMap(classStr);
                end
            end
        end
    end

    methods (Static)
        function b = ShouldIgnoreLabels(labels)
            b = (labels == Labeling.kIgnoreClass | labels == Labeling.kInvalidClass);
        end

        function classStr = StringForClass(class, classNames)
            if class == Labeling.kNullClass
                classStr = Labeling.kNullClassStr;
            else
                classStr = classNames{class};
            end
        end
    end

    methods (Access = private)
        function createLabeling(obj,classNames)
            if ~isempty(classNames)
                nClasses = length(classNames);
                obj.classesMap = containers.Map(classNames,int8(1:nClasses));
                obj.classesMap(Labeling.kIgnoreStr) = Labeling.kIgnoreClass;
                obj.classesMap(Labeling.kNullClassStr) = Labeling.kNullClass;
            end
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
classdef CI_HS < PROBLEM
% <single> <real> <large/none> <multitask>
% Multitasking problem (Griewank function + Rastrigin function)
% SubD --- 50,50 --- Number of decision variables of each task

%------------------------------- Reference --------------------------------
% K. K. Bali, Y. Ong, A. Gupta, and P. S. Tan, Multifactorial evolutionary
% algorithm with online transfer parameter estimation: MFEA-II, IEEE
% Transactions on Evolutionary Computation, 2020, 24(1): 69-83.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        SubD;   % Number of decision variables of each task
        L1;   	% Low bounds of the first task
        L2;   	% Low bounds of the second task
        U1;   	% Upper bounds of the first task
        U2;   	% Upper bounds of the second task
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.SubD     = obj.ParameterSet([50,50]);
            obj.M        = 1;
            obj.D        = max(obj.SubD) + 1;
            obj.L1       = zeros(1,obj.SubD(1)) - 600;
            obj.U1       = zeros(1,obj.SubD(1)) + 600;
            obj.L2       = zeros(1,obj.SubD(2)) - 5.12;
            obj.U2       = zeros(1,obj.SubD(2)) + 5.12;
            obj.lower    = [zeros(1,obj.D-1),1];
            obj.upper    = [ones(1,obj.D-1),length(obj.SubD)];
            obj.encoding = [ones(1,obj.D-1),2];
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            PopObj = zeros(size(PopDec,1),1);
            for i = 1 : size(PopDec,1)
                if PopDec(i,end) == 1       % Task 1
                    x1 = obj.L1 + PopDec(i,1:obj.SubD(1)).*(obj.U1-obj.L1);
                    PopObj(i) = 1/4000*sum(x1.^2,2) - prod(cos(x1./sqrt(repmat(1:size(x1,2),size(x1,1),1))),2) + 1;
                elseif PopDec(i,end) == 2   % Task 2
                    x2 = obj.L2 + PopDec(i,1:obj.SubD(2)).*(obj.U2-obj.L2);
                    PopObj(i) = sum(x2.^2-10*cos(2*pi*x2)+10,2);
                end
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% Weak-Texture Mask generation
%
% msk = WeakTextureMask(img, patchsize, th)
%
%Output parameters
% msk: weak-texture mask. 0 and 1 represent non-weak-texture and weak-texture regions, respectively
%
%
%Input parameters
% img: input single image
% th: threshold which is output of NoiseLevel
% patchsize (optional): patch size (default: 7)
%
%Example:
% img = double(imread('img.png'));
% patchsize = 7;
% [nlevel th] = NoiseLevel(img, patchsize);
% msk = WeakTextureMask(img, patchsize, th);
% imwrite(uint8(msk*255), 'msk.png');
% version: 20150203

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Noise Level Estimation:                                       %
%                                                               %
% Copyright (C) 2012-2015 Masayuki Tanaka. All rights reserved. %
%                    mtanaka@ctrl.titech.ac.jp                  %
%                                                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function msk = WeakTextureMask(img, th, patchsize)

if( ~exist('patchsize', 'var') )
    patchsize = 7;
end

kh = [-1/2,0,1/2];
imgh = imfilter(img,kh,'replicate');
imgh = imgh(:,2:size(imgh,2)-1,:);
imgh = imgh .* imgh;

kv = kh';
imgv = imfilter(img,kv,'replicate');
imgv = imgv(2:size(imgv,1)-1,:,:);
imgv = imgv .* imgv;

numRows     = size(img, 1);
numCols     = size(img, 2);
numChannels = size(img, 3);

msk = zeros(numRows, numCols, numChannels);

for cha=1:numChannels
	m = im2col(img(:,:,cha),[patchsize patchsize]);
	m = zeros(size(m));
	Xh = im2col(imgh(:,:,cha),[patchsize patchsize-2]);
	Xv = im2col(imgv(:,:,cha),[patchsize-2 patchsize]);

	Xtr = sum(vertcat(Xh,Xv));

	p = (Xtr<th(cha));
	ind = 1;
	for col=1:numCols-patchsize+1
		for row=1:numRows-patchsize+1
			if( p(ind) > 0 )
				msk(row:row+patchsize-1, col:col+patchsize-1, cha) = 1;
			end
			ind = ind + 1;
		end
	end

end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Example code for estimating the HRF using the Inverse-Logit Model, a
% Finte Impluse Response Model and the Canonical HRF with 2 derivatives.
% Also the code illustrates our code for detecting model misspecification.
%
% By Martin Lindquist and Tor Wager
% Created  10/02/09
% Last edited 05/20/13
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Load time course
%

mypath = which('ilogit');
if isempty(mypath), error('Cannot find directory with ilogit.m and other functions. Not on path?'); end
[mydir] = fileparts(mypath)

load(fullfile(mydir,'timecourse'))

tc = (tc- mean(tc))/std(tc);
len = length(tc);

%% Or: create your own
[xBF] = spm_get_bf(struct('dt', .5, 'name', 'hrf (with time and dispersion derivatives)', 'length', 32));
clear Xtrue
for i = 1:1, xx = conv(xBF.bf(:,i), [1 1 1 1 1 1 ]');
    Xtrue(:, i) = xx(1:66);
end
for i = 2:3, xx = conv(xBF.bf(:,i), [1 1]');
    Xtrue(:, i) = xx(1:66);
end
hrf = Xtrue * [1 .3 .2]';
xsecs = 0:.5:32;

hrf = [ 0; 0; hrf];
hrf = hrf(1:length(xsecs));
hrf = hrf ./ max(hrf);
figure; plot(xsecs, hrf, 'k')
%hrf = hrf(1:4:end); % downsample to TR, if TR is > 0.5

R = randperm(640); R = sort(R(1:36));
Run = zeros(640,1);
for i=1:length(R), Run(R(i)) = 1; end;
true_sig = conv(Run, hrf);
true_sig = true_sig(1:640);

% tc_noise = noise_arp(640, [.7 .2]);
% tc = true_sig + 0.1 * tc_noise;
tc = true_sig;
%figure; plot(tc);

Runc{1} = Run;

%%

create_figure; subplot(3,1,1); han = plot(tc);
title('Sample time course'); drawnow

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Settings
%

TR = 0.5;
%T = round(30/TR);
T = 30;
t = 1:TR:T;                        % samples at which to get Logit HRF Estimate
FWHM = 4;                       % FWHM for residual scan
pval = 0.01;
df = 600;
alpha = 0.001;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Create stick function (sample event onsets)
% Variable R contains onset times
% Variable Run contains stick (a.k.a. delta or indicator) function

% R = [3, 21, 56, 65, 109, 126, 163, 171, 216, 232, 269, 282, 323, 341, 376, 385, 429, 446, 483, 491, 536, 552, 589, 602];
% Run = zeros(640,1);
% for i=1:length(R), Run(R(i)) = 1; end;
%

try
    hold on;
    hh = plot_onsets(R,'k',-3,1, 1);
    drawnow
catch
    disp('Couldn''t find function to add onset sticks to plot. Skipping.')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fit HRF using IL-function

% Choose mode (deterministic/stochastic)

mode = 0;   % 0 - deterministic aproach
            % 1 - simulated annealing approach
            % Please note that when using simulated annealing approach you
            % may need to perform some tuning before use.

[h1, fit1, e1, param] = Fit_Logit2(tc,TR,Runc,T,mode);
[pv sres sres_ns1] = ResidScan(e1, FWHM);
[PowLoss1] = PowerLoss(e1, fit1, (len-7) , tc, TR, Runc, alpha);

hold on; han(2) = plot(fit1,'r');

disp('Summary: IL_function');

disp('Amplitude:'); disp(param(1));
disp('Time-to-peak:'); disp(param(2)*TR);
disp('Width:'); disp(param(3)*TR);

disp('MSE:'); disp((1/(len-1)*sum(e1.^2)));
disp('Mis-modeling:'); disp(pv);
disp('Power Loss:'); disp(PowLoss1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fit HRF using FIR-model

% Choose mode (FIR/sFIR)

mode = 1;   % 0 - FIR
            % 1 - smooth FIR

[h2, fit2, e2, param] = Fit_sFIR(tc,TR,Runc,T,mode);
[pv sres sres_ns2] = ResidScan(e2, FWHM);
[PowLoss2] = PowerLoss(e2, fit2, (len-T) , tc, TR, Runc, alpha);

hold on; han(3) = plot(fit2,'g');

disp('Summary: FIR');

disp('Amplitude'); disp(param(1));
disp('Time-to-peak'); disp(param(2)*TR);
disp('Width'); disp(param(3)*TR);

disp('MSE:'); disp((1/(len-1)*sum(e2.^2)));
disp('Mis-modeling'); disp(pv);
disp('Power Loss:'); disp(PowLoss2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fit HRF using Canonical HRF + 2 derivatives

p=1;

[h3, fit3, e3, param, info] = Fit_Canonical_HRF(tc,TR,Runc,30,p);
[pv sres sres_ns3] = ResidScan(e3, FWHM);
[PowLoss3] = PowerLoss(e3, fit3, (len-p) , tc, TR, Runc, alpha);

hold on; han(4) = plot(fit3,'m');

legend(han,{'Data' 'IL' 'sFIR' 'DD'})

disp('Summary: Canonical + 2 derivatives');

disp('Amplitude'); disp(param(1));
disp('Time-to-peak'); disp(param(2)*TR);
disp('Width'); disp(param(3)*TR);

disp('MSE:'); disp((1/(len-1)*sum(e3.^2)));
disp('Mis-modeling'); disp(pv);
disp('Power Loss:'); disp(PowLoss3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%figure;
%%

subplot(3,2,5); hold on;
plot(xsecs, hrf, 'k')
xsecs1 = xsecs(1:length(h1));
han2 = plot(xsecs1, h1,'r');
xsecs2 = xsecs(1:length(h2));
han2(2) = plot(xsecs2, h2,'g');
xsecs3 = xsecs(1:length(h3));
han2(3) = plot(xsecs3, h3,'m');
legend(han2,{'IL' 'sFIR' 'DD'})
title('Estimated HRF');

subplot(3,1,2); hold on;
hh = plot_onsets(R,'k',-3,1);
drawnow

han3 = plot(sres_ns1,'r');
hold on; han3(2) = plot(sres_ns2,'g');
hold on; han3(3) = plot(sres_ns3,'m');
hold on; plot((1:len),zeros(len,1),'--k');
legend(han3,{'IL' 'sFIR' 'DD'})
title('Mis-modeling (time course)');

subplot(3,2,6); hold on;

[s1] = Fit_sFIR(sres_ns1,TR,Runc,T,0);
[s2] = Fit_sFIR(sres_ns2,TR,Runc,T,0);
[s3] = Fit_sFIR(sres_ns3,TR,Runc,T,0);

han4 = plot(s1(1:T),'r');
hold on; han4(2) = plot(s2(1:T),'g');
hold on; han4(3) = plot(s3(1:T),'m');
hold on; plot((1:T),zeros(T,1),'--k');
legend(han4,{'IL' 'sFIR' 'DD'})
title('Mis-modeling (HRF)');

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function K = rbfhKernCompute(kern, x1, x2)

% RBFHKERNCOMPUTE Compute the RBFH kernel given the parameters and X.
% FORMAT
% DESC computes the kernel parameters for the radial basis function heat
% kernel given inputs associated with rows and columns.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x1 : the input matrix associated with the rows of the kernel.
% ARG x2 : the input matrix associated with the columns of the kernel.
% RETURN K : the kernel matrix computed at the given points.
%
% FORMAT
% DESC computes the kernel matrix for the radial basis function heat
% kernel given a design matrix of inputs.
% ARG kern : the kernel structure for which the matrix is computed.
% ARG x1 : input data matrix in the form of a design matrix.
% RETURN K : the kernel matrix computed at the given points.
%
% SEEALSO : rbfhKernParamInit, kernCompute,
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if nargin < 3
    x2 = x1;
end

if size(x1, 2) ~= 2 || size(x2, 2) ~= 2
    error('Input can only have two columns');
end

% Split the domain into time domain and spatial domain and account for
% missing values. If there are no missing values the computation of the
% kernel is a pointwise prodruct, otherwise it is a kronecker product.
t1 = x1(x1(:,1)~=Inf,1);
t2 = x2(x2(:,1)~=Inf,1);
s1 = x1(x1(:,2)~=Inf,2);
s2 = x2(x2(:,2)~=Inf,2);
if (length(t1) == length(s1)) && (length(t2) == length(s2))
    ut1 = unique(t1);
    ut2 = unique(t2);
    us1 = unique(s1);
    us2 = unique(s2);
    if (length(ut1)*length(us1) == length(t1)) && ...
        (length(ut2)*length(us2) == length(t2))
        t1 = ut1; s1 = us1; t2 = ut2; s2 = us2;
        isPointwise = false;
    else
        isPointwise = true;
    end
else
    isPointwise = false;
end

kern.rbf.inverseWidth = kern.inverseWidthTime;
Kt = rbfKernCompute(kern.rbf, t1, t2);
kern.rbf.inverseWidth = kern.inverseWidthSpace;
Ks = rbfKernCompute(kern.rbf, s1, s2);

if isPointwise
   K = Kt.*Ks;
else
   K = kron(Kt, Ks);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
close all;
clearvars;
clc;
rng default;
N = 32;
K = 4;
gen = spx.data.synthetic.SparseSignalGenerator(N, K);
rep =  gen.uniform();
figure;
stem(rep, '.');
export_fig images/demo_sparse_uniform_1.png -r120 -nocrop;

author: indigits

repo: sparse-plex


text:
function vl_figaspect(a)
% VL_FIGASPECT  Set figure aspect ratio
%   VL_FIGASPECT(A) sets the current figure aspect ratio to A. It
%   leaves the length of the shortest side unaltered. Both the
%   screen and paper positions are affected.
%
%   See also: VL_PRINTSIZE(), VL_HELP().

% Authors: Andrea Vedaldi

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

pos = get(gcf, 'Position') ;

if a >= 1
  pos(3) = a * pos(4) ;
else
  pos(4) = pos(3) / a ;
end

set(gcf,'Position',pos) ;

pos = get(gcf, 'PaperPosition') ;

if a >= 1
  pos(3) = a * pos(4) ;
else
  pos(4) = pos(3) / a ;
end

set(gcf,'PaperPosition', pos) ;

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function model = robTwoDynamicsSetLatentValues(model, X)

% ROBTWODYNAMICSSETLATENTVALUES Set the latent values inside the model.

% FGPLVM

model = robOneDynamicsSetLatentValues(model, X);


text:
function [offsets,WI,parents,I] = bone_forest(C,E,P,BE)
  % [offsets,WI,parents] = bone_forest(C,E,P,BE)
  %
  % See also: writeBF.m

  if nargin == 2 || isempty(BE)
    BE = E;
    P = [];
  end
  [BE,I] = sortrows(BE);
  % point properties
  P_parents = zeros(numel(P),1);
  P_offsets = C(P,:);
  P_WI = (1:numel(P))';
  % bone properties
  % bones form tree(s)
  assert(numel(unique(BE(:,2))) == numel(BE(:,2)));
  % 0 if not dest otherwise index in list
  index_map = zeros(max(BE(:)),1);
  index_map(BE(:,2)) = 1:size(BE,1);
  roots = setdiff(BE(:,1),BE(:,2));
  % add roots to index map
  index_map(roots) = size(BE,1)+(1:numel(roots));
  % parents after adding roots
  BE_parents = index_map(BE(:,1));
  % roots have no mother or father
  BE_parents = [BE_parents; zeros(numel(roots),1)];
  BE_WI = [1:size(BE,1) zeros(1,numel(roots))]';
  BE_offsets = bone_offsets(C([BE(:,2); roots],:),BE_parents);
  % we'll first put "point" handles
  % then bone's (recall that in BF each joint needs an entry, including roots)
  parents = [P_parents;numel(P).*(BE_WI~=0)+BE_parents];
  WI = [P_WI;numel(P).*(BE_WI~=0)+BE_WI];
  offsets = [P_offsets;BE_offsets];
end


text:
function selective_average_interactive_view_init(imgs, onsets, scans_per_sess, TR, hp_length, t, basepts, plotstes)
% :Usage:
% ::
%
%    selective_average_interactive_view_init(imgs, onsets, 20, 1:2, 1);
%
% Initialize point-and-click data extraction and plotting of selective
% averages
%
% :Inputs:
%
%   **imgs:**
%        cell array of image files, one cell per subject
%
%        imgs is cell array of time series image names for each subject
%
%   **onsets:**
%        cell array of onsets
%
%        Each cell contains its own cell array, with one cell per
%        condition.
%
%   **t:**
%        time points in average to estimate
%
%   **basepts:**
%        indices of baseline points to subtract from individual averages
%
%   **plotstes:**
%        plot standard errors: 1 or 0

    disp('Initializing Selective Average interactive viewer') % need to do the stuff below only once

    disp('-------------------------------------------------')

    %     disp('Loading mediation_SETUP.mat')
    %     [SETUP, imgs, wh_is_image, name] = load_mediation_setup;

    N = length(imgs);

    % Setting up filtering
    % ------------------------------------------------------
    fprintf('Setting up high-pass filtering (dummy params for first 2 time points per run): ');
    % %     I = cell(1, N);
    % %     S = cell(1, N);
    tic

    for i = 1:N
        [tmp, I, S] = hpfilter(rand(sum(scans_per_sess), 1), TR, hp_length, scans_per_sess, [], 1:2);
    end

    fprintf(' %3.0f s\n', toc)

    % ------------------------------------------------------
    V = cell(1, N);
    fprintf('Mapping image volumes to memory: ');
    tic
    for i = 1:N
        fprintf('%3.0f ', i);
        V{i} = spm_vol(imgs{i});
    end

    fprintf(' %3.0f s\n', toc)

    % ------------------------------------------------------

    disp('Registering graphics callback with SPM registry')

    callback_handle = @(str, pos, reg, hReg) selavg_interactive_callback_wrapper(str, pos, reg, hReg);

    hSpmFig = spm_figure('GetWin', 'Graphics');

    hReg = uicontrol(hSpmFig, 'Style', 'Text', 'String', 'InteractiveViewer hReg', ...
        'Position', [100 200 100 025], 'Visible', 'Off', ...
        'FontName', 'Times', 'FontSize', 14, 'FontWeight', 'Bold', ...
        'HorizontalAlignment', 'Center');

    hReg = spm_XYZreg('InitReg', hReg, V{1}(1).mat, V{1}(1).dim(1:3)');
    spm_XYZreg('Add2Reg', hReg, 0, callback_handle);
    spm_orthviews('Register', hReg);

    %fh = findobj('Tag', 'Graphics');
    %set(fh, 'WindowButtonUpFcn', callback_handle);

    disp('Ready!')

    % inline

    function selavg_interactive_callback_wrapper(str, pos, reg, hReg)

        switch str
            case 'SetCoords'
                selavg_interactive_callback(pos, onsets, V, t, basepts, plotstes, S, scans_per_sess, I);

            otherwise
                disp('Unknown callback command from spm_XYZreg');
        end

    end

end

function selavg_interactive_callback(pos, onsets, V, t, basepts, plotstes, S, scans_per_sess, I)
    %% Done each time you click

    % %     %pos = spm_orthviews('Pos');
    % %     vox = (mm2voxel(pos', V{1}(1)))';
    % %     vox(4) = 1;
    % %
    % %     fprintf('Position: %3.0f %3.0f %3.0f\n', pos(1), pos(2), pos(3));
    % %     go = input('Plot selective averages for this voxel? (1 or 0) : ');

    % %     if go
    % % %         % Load data for this voxel
    % % %         % -----------------------------------------------------
    % % %         N = length(V);
    % % %
    % % %         fprintf('Loading voxel data for all subjects: ');
    % % %         for i = 1:N, fprintf('%3.0f ', i); braindata{i} = spm_get_data(V{i}, vox); end

    % Load data and get selective averages for each subject
    % -----------------------------------------------------
    [group_avgs, group_stes, subject_avgs, subject_stes, braindata] = selective_average_group( ...
        V, onsets, pos, 't', t, 'basepts', basepts, 'plotstes', plotstes, 'S', S, 'scans', scans_per_sess, 'I', I);

    if ~isempty(group_avgs)

        disp('-------------------------------------------------------');
        disp('Assigning group_avgs and stes in base workspace.');
        disp('(and subject_avgs and subject_stes)');
        disp('Assigning braindata to base workspace: Timeseries for each subject');
        disp('-------------------------------------------------------');
        disp(' ')

        assignin('base', 'group_avgs', group_avgs);
        assignin('base', 'group_stes', group_stes);
        assignin('base', 'subject_avgs', subject_avgs);
        assignin('base', 'subject_stes', subject_stes);
        assignin('base', 'braindata', braindata);

    end
    % %     end
end

function [SETUP, imgs, wh_is_image, name] = load_mediation_setup

    SETUP = [];
    imgs = [];

    fname = [pwd filesep 'mediation_SETUP.mat'];
    if exist(fname,'file')
        load(fname);

        % try to find names (single level)
        if exist('SETUP','var') && isfield(SETUP, 'M') && ischar(SETUP.M)
            imgs = SETUP.M;
            name = 'From mediation_SETUP SETUP.M';
            wh_is_image = 'M';

        elseif exist('SETUP','var') && isfield(SETUP, 'X') && ischar(SETUP.X)
            imgs = SETUP.X;
            name = 'From mediation_SETUP SETUP.X';
            wh_is_image = 'X';

        elseif exist('SETUP','var') && isfield(SETUP, 'Y') && ischar(SETUP.Y)
            imgs = SETUP.Y;
            name = 'From mediation_SETUP SETUP.Y';
            wh_is_image = 'Y';

        end

        % try to find names: multi-level
        if isfield(SETUP, 'data')
            switch SETUP.cmdstring
                case 'Search for mediators'
                    imgs = SETUP.data.M;
                    name = 'Multilevel, from SETUP.data.M';
                    wh_is_image = 'M';

                case 'Search for indirect influences'
                    imgs = SETUP.data.X;
                    name = 'Multilevel, from SETUP.data.X';
                    wh_is_image = 'X';

                case 'Search for mediated outcomes'
                    imgs = SETUP.data.Y;
                    name = 'Multilevel, from SETUP.data.Y';
                    wh_is_image = 'Y';

                otherwise
                    error('Unknown cmdstring: "%s".', cmdstring);
            end

            if ~iscell(imgs) || ~ischar(imgs{1})
                imgs = []; % invalid data here
            end

        end

        if isempty(imgs)
            fprintf(1,'Could not find image list.\n');
        end

    else
        fprintf(1,'Go to valid mediation directory with SETUP.mat to use interactive plotting.\n');
    end

end


text:
function cellInfo = Plant_States_Bus(varargin)
% PLANT_STATES_BUS returns a cell array containing bus object information
%
% Optional Input: 'false' will suppress a call to Simulink.Bus.cellToObject
%                 when the MATLAB file is executed.
% The order of bus element attributes is as follows:
%   ElementName, Dimensions, DataType, SampleTime, Complexity, SamplingMode, DimensionsMode, Min, Max, DocUnits, Description

suppressObject = false;
if nargin == 1 && islogical(varargin{1}) && varargin{1} == false
    suppressObject = true;
elseif nargin > 1
    error('Invalid input argument(s) encountered');
end

cellInfo = { ...
  { ...
    'Plant_States_Bus', ...
    '', ...
    '', ...
    'Auto', ...
    '-1', {...
{'timestamp', 1, 'uint32', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('ms'), ''}; ...
{'phi', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('rad'), ''}; ...
{'theta', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('rad'), ''}; ...
{'psi', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('rad'), ''}; ...
{'rot_x_B', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('rad/s'), ''}; ...
{'rot_y_B', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('rad/s'), ''}; ...
{'rot_z_B', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('rad/s'), ''}; ...
{'acc_x_O', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('m/s2'), ''}; ...
{'acc_y_O', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('m/s2'), ''}; ...
{'acc_z_O', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('m/s2'), ''}; ...
{'vel_x_O', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('m/s'), ''}; ...
{'vel_y_O', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('m/s'), ''}; ...
{'vel_z_O', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('m/s'), ''}; ...
{'x_R', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('m'), ''}; ...
{'y_R', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('m'), ''}; ...
{'h_R', 1, 'single', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('m'), ''}; ...
{'lat', 1, 'double', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('rad'), ''}; ...
{'lon', 1, 'double', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('rad'), ''}; ...
{'alt', 1, 'double', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('m'), ''}; ...
{'lat_0', 1, 'double', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('rad'), ''}; ...
{'lon_0', 1, 'double', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('rad'), ''}; ...
{'alt_0', 1, 'double', -1, 'real', 'Sample', 'Fixed', [], [], sprintf('m'), ''}; ...
    } ...
  } ...
}';

if ~suppressObject
    % Create bus objects in the MATLAB base workspace
    Simulink.Bus.cellToObject(cellInfo)
end

author: Firmament-Autopilot

repo: FMT-Model


text:
function [favar]=favar_irf(favar,IRFperiods,endo,IRFt,strctident,pref)

% we arranged everything before including transformation
    % pick the relevant shocks (plotXshock)
    favar_irf_estimates=favar.irf_estimates(:,favar.IRF.plotXshock_index);

    % save in favar structure, ready to plot
    favar.IRF.favar_irf_estimates=favar_irf_estimates;

    % plot the IRFs of the factors
    if favar.IRF.plot==1
    irf_favar=figure;
    set(irf_favar,'name','approximate impulse response functions (FAVAR)');
    numIRFs=favar.npltX*favar.IRF.npltXshck; %total number of IRFs to print
    numcol=5;
    numrow=ceil(numIRFs/numcol); % maximum number of rows

    for ii=1:numIRFs
        subplot(numrow,numcol,ii);
        hold on
        Xpatch=[(1:IRFperiods) (IRFperiods:-1:1)];
        Ypatch=[favar.IRF.favar_irf_estimates{ii}(1,:) fliplr(favar.IRF.favar_irf_estimates{ii}(3,:))];
        IRFpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
        set(IRFpatch,'facealpha',0.5);
        set(IRFpatch,'edgecolor','none');
        plot(favar.IRF.favar_irf_estimates{ii}(2,:),'Color',[0.4 0.4 1],'LineWidth',2);
        plot([1,IRFperiods],[0 0],'k--');
        hold off
%         minband=min(favar.IRF.irf_factors_final_plotX{ii}(1,:));
%         maxband=max(favar.IRF.irf_factors_final_plotX{ii}(3,:));
%         space=maxband-minband;
%         Ymin=minband-0.2*space;
%         Ymax=maxband+0.2*space;
        set(gca,'XLim',[1 IRFperiods],'FontName','Times New Roman'); %,'YLim',[Ymin Ymax]
        % top labels
   if ii<=favar.IRF.npltXshck
       if IRFt==1||IRFt==2||IRFt==3
       printlabels=endo{favar.IRF.plotXshock_index(1,ii),1};
       elseif IRFt==4||IRFt==5||IRFt==6
       printlabels=strctident.signreslabels{strctident.signreslabels_shocksindex(ii,1)};
        end
   title(printlabels,'FontWeight','normal','interpreter','none');
   end
% side labels
   if rem((ii-1)/favar.IRF.npltXshck,1)==0
   ylabel(favar.informationvariablestrings{1,favar.plotX_index((ii-1)/favar.IRF.npltXshck+1)},'FontWeight','normal','interpreter','none');
   end
    end
    % top supertitle
    ax=axes('Units','Normal','Position',[.11 .075 .85 .88],'Visible','off');
    set(get(ax,'Title'),'Visible','on')
    title('Shock:','FontSize',11,'FontName','Times New Roman','FontWeight','normal');
    % side supertitle
    ylabel('approximate response of:','FontSize',12,'FontName','Times New Roman','FontWeight','normal');
    set(get(ax,'Ylabel'),'Visible','on')
    end

    %% save in excel
    excelrecord4favar

author: european-central-bank

repo: BEAR-toolbox


text:
function i4_gcd_test ( )

%*****************************************************************************80
%
%% I4_GCD_TEST tests I4_GCD.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  test_num = 7;

  i_test = [ 36, 49, 0, 12, 36, 1, 91 ];
  j_test = [ 30, -7, 71, 12, 49, 42, 28 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4_GCD_TEST\n' );
  fprintf ( 1, '  I4_GCD computes the greatest common factor\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I     J   I4_GCD\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num
    i = i_test(test);
    j = j_test(test);
    fprintf ( 1, '  %6d  %6d  %6d\n', i, j, i4_gcd ( i, j ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [label, R] = mixGaussPred(model, X)
% Predict label and responsibility for Gaussian mixture model.
% Input:
%   X: d x n data matrix
%   model: trained model structure outputed by the EM algirthm
% Output:
%   label: 1 x n cluster label
%   R: k x n responsibility
% Written by Mo Chen (sth4nth@gmail.com).
mu = model.mu;
Sigma = model.Sigma;
w = model.w;

n = size(X,2);
k = size(mu,2);
logRho = zeros(n,k);

for i = 1:k
    logRho(:,i) = loggausspdf(X,mu(:,i),Sigma(:,:,i));
end
logRho = bsxfun(@plus,logRho,log(w));
T = logsumexp(logRho,2);
logR = bsxfun(@minus,logRho,T);
R = exp(logR);
[~,label(1,:)] = max(R,[],2);

function y = loggausspdf(X, mu, Sigma)
d = size(X,1);
X = bsxfun(@minus,X,mu);
[U,p]= chol(Sigma);
if p ~= 0
    error('ERROR: Sigma is not PD.');
end
Q = U'\X;
q = dot(Q,Q,1);  % quadratic term (M distance)
c = d*log(2*pi)+2*sum(log(diag(U)));   % normalization constant
y = -(c+q)/2;

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function parsedGPR = GPRparser(model, getCNFSets)
% Maps the GPR rules of the model to a specified format that is used by
% the model extraction methods
%
% USAGE:
%   parsedGPR = GPRparser(model)
%
% INPUT:
%   model:       cobra model structure
%
% OPTIONAL INPUT:
%    getCNFSets:    whether to get the CNF sets (true) or DNF sets (false).
%                   DNF sets represent functional enzyme complexes, while
%                   CNF sets represent the possible subunits of a complex.
%                   (default: false , i.e. DNF sets)
%
% OUTPUT:
%   parsedGPR:   cell matrix containing parsed GPR rule
%
% AUTHORS: Thomas Pfau & Anne Richelle, May 2017

if ~exist('getCNFSets','var')
    getCNFSets = false;
end

parsedGPR = {};
fp = FormulaParser();
for i = 1:numel(model.rxns)
    if ~isempty(model.rules{i})
        head = fp.parseFormula(model.rules{i});
        currentSets = head.getFunctionalGeneSets(model.genes,getCNFSets)';
        parsedGPR{i}=currentSets;
    else
        parsedGPR{i}={''};
    end
end
parsedGPR=parsedGPR';
end

author: opencobra

repo: cobratoolbox


text:
%%% void show_3dvideo(3dmatrix)
%
function show_3dvideo(V)
  for i = 1 : size(V,3)
    frame = V(:,:,i);
    imshow(frame,[],'InitialMagnification','fit');
    disp(i);
    pause(0.01);
  end
end

author: andrewssobral

repo: lrslibrary


text:
% Demo script for software from the paper
% 'Generating spike-trains with specified correlations',
%  Macke et al., submitted for publication, 2008
%
% www.kyb.mpg.de/bethgegroup/code/efficientsampling
%
% Instructions:
%  - Change to the code directory and run demo.m
%  - The script will take you through the functions .
%  - After each step, the script will pause.
%  - To continue, hit any button.
%  - Read along in the demo.m file to follow what's happening.
%
% If you have questions, feel free to email us.

f = 1;
fprintf('\nDemo program for GENERATING SPIKE TRAINS WITH SPECIFIED CORRELATIONS, Macke et al.\n')

%% Step 1: Generating correlated binary variables (2D example)
mu = [.4,.3]';      % set the mean
v = mu.*(1-mu);     % calculate the variances
C = diag(v);        % build covariance matrix
C(1,2) = .1;
C(2,1) = .1;

[S,g,L] = sampleDichGauss01(mu,C,1e5);   % generate samples from the DG model

muHat = mean(S,2);  % estimate mean
CHat = cov(S');     % estimate covariance

fprintf('\nStep 1: Generating correlated binary variables (section 2.1)\n\n')
fprintf('Target mean X1:  %.2f     Estimated mean X1:  %.3f\n',mu(1),muHat(1))
fprintf('Target mean X2:  %.2f     Estimated mean X2:  %.3f\n',mu(2),muHat(2))
fprintf('Target cov X1X2: %.2f     Estimated cov X1X2: %.3f\n',C(1,2),CHat(1,2))

disp('To proceed to compare histograms, hit any key...')
pause

%% Step 2: How much do correlations change the distribution when means are
% assumed to be equal?

fprintf('\nStep 2: Effect of correlations\n\n')
fprintf('Computing...\n')

% First, we generate a mean vector in ten dimensions. We assume that all
% means are approximately .5, but we add some Gaussian noise to the network
% is not to uniform. We set the covariances so that the resulting
% correlations are sizable, but not to large (~0.1).

mu = ones(10,1) * .2 + randn(10,1)*0.05;
v = mu.*(1-mu);
C = ones(10,10) * .02;
C(eye(size(C))==1) = v;

% We find the histogram of the independent distribution with this mean
% vector P(x) = PROD_i P(X_i=1) and the DG distribution with the same mean
% and the covariances as defined above.

h1 = binHistIndep(mu);    % returns the histogram of an independet distribution
S = sampleDichGauss01(mu,C,1e5);
h2 = binHist(S);          % creates a histogram of the binary patterns in S
h2 = h2 / sum(h2);

% When we compare them, we see that some patterns occur much more often
% in the correlated then in the independent distribution. This a clear
% indication of the strength of the correlations. The "stripes" coincide
% with patterns with a different number of active neurons (starting at 0 in
% the upper left corner).

figure(f)
loglog(h1,h2,'k.')
xlabel('P(X) in independent model')
ylabel('P(X) in DG model')
title('Step 2: What effect do correlations have?')
axis square

f = f+1;

disp('To proceed to generate Poisson RVs, hit any key...')
pause

%% Step 3: Correlated Poisson from a binary process (3.1)
% First, we generate samples from a correlated Poisson distribution with
% specified mean and covariance matrix. The covariance between the two
% distributions is positive and quite strong. The construction follows the
% description in section 3.1, i.e. we first find a binary process with the
% correct parameters and sum over it. This means we create
% a joint Poisson distribution.

fprintf('\nStep 3: Correlated Poisson from a binary process (3.1)\n\n')
fprintf('Computing...\n')

mu = [7,9]';      % set the mean
C = [7 3;3 9];    % set an admissable covariance matrix

S = sampleCovPoisson(mu,C,10000);   % generate sample via discretized Gaussian

% We obtain a set of samples S. Next, we verify that the samples have
% desired mean and covariance structure. We see that this is indeed the case.

muHat = mean(S,2);
CHat = cov(S');

fprintf('Target mean X1:  %.2f     Estimated mean X1:  %.3f\n',mu(1),muHat(1))
fprintf('Target mean X2:  %.2f     Estimated mean X2:  %.3f\n',mu(2),muHat(2))
fprintf('Target cov X1X2: %.2f     Estimated cov X1X2: %.3f\n',C(1,2),CHat(1,2))

disp('To plot the marginals and joint histogram, hit any key...')
pause
fprintf('Computing...\n')

figure (f)

% Now, we look at the marginal distribution we obtain from our sampling
% procedure. For comparison, we also plot the true 1D Poisson distributions
% as specified by the mean of the marginal (dotted). We see that they are
% very close to one another.

h1 = histc(S(1,:),0:30); h1 = h1/sum(h1);
h2 = histc(S(2,:),0:30); h2 = h2/sum(h2);

subplot(1,2,1)
plot(0:length(h1)-1,h1,'-r','markersize',5), hold on
plot(0:length(h1),poisspdf(0:length(h1),mu(1)),'.r','markersize',5)
plot(0:length(h2)-1,h2,'-g','markersize',5)
plot(0:length(h1),poisspdf(0:length(h1),mu(2)),'.g','markersize',5)
t = legend('X_1','X_1 Poisson','X_2','X_2 Poisson');
set(t,'box','off')
axis square
xlabel('X'), ylabel('P(X)')
title(sprintf('Step 3: Correlated Poisson, Positive Correlation\nMarginal Distributions with Poisson\n distribution with correct mean'))

% Finally, we convince ourselves that the two variables are indeed
% positively correlated by looking at the 2D joint histogram. We can see
% that although the marginals are perfect Poissonian, most samples are
% concentrated along the main diagonal, indicating the positive
% correlation.

subplot(1,2,2)
hh = EstimateDiscreteJoint(S);
imagesc(hh)
axis square
title('')
xlabel('X_1'), ylabel('X_2')

disp('To proceed to positively correlated Poisson variables with another method, hit any key...')
pause

f = f+1;

%% Step 4: Generating correlated Poisson variables via the method
% in section 3.3 for positive correlations
%
% Again, we generate samples from a correlated Poisson distribution with
% specified mean and covariance matrix. The covariance between the two
% distributions is positive and quite strong. The construction follows now
% section 3.3, so we truncate a Gaussian to obtain a distribution with
% Poisson marginals and a given covariance.

fprintf('\nStep 4: Generating positively correlated Poisson variables (section 3.3)\n\n')
fprintf('Computing...\n')

mu = [7,9]';      % set the mean
C = [7 3;3 9];    % set an admissable covariance matrix

[S,gamma,Lambda,joints2D] = DGPoisson(mu,C,1e5);   % generate sample via discretized Gaussian

% We obtain a set of samples S, the parameters of the hidden Gaussian
% variable (gamma and Lambda) and a structure containing the marginal
% distributions as well as the 2D joint distribution

% Next, we verify that the samples have desired mean and covariance
% structure. We see that this is indeed the case.

muHat = mean(S,2);
CHat = cov(S');

fprintf('Target mean X1:  %.2f     Estimated mean X1:  %.3f\n',mu(1),muHat(1))
fprintf('Target mean X2:  %.2f     Estimated mean X2:  %.3f\n',mu(2),muHat(2))
fprintf('Target cov X1X2: %.2f     Estimated cov X1X2: %.3f\n',C(1,2),CHat(1,2))

disp('To plot the marginals and joint histogram, hit any key...')
pause

figure(f)

% Now, we look at the marginal distribution we obtain from our sampling
% procedure. For comparison, we also plot the true 1D Poisson distributions
% as specified by the mean of the marginal (dotted). We see that they are
% very close to one another.

subplot(2,2,1)
h1 = joints2D{1,1};   % marginal of X_1
h2 = joints2D{2,2};   % marginal of X_2
plot(0:length(h1)-1,h1,'-r','markersize',5), hold on
plot(0:length(h1),poisspdf(0:length(h1),mu(1)),'.r','markersize',5)
plot(0:length(h2)-1,h2,'-g','markersize',5)
plot(0:length(h1),poisspdf(0:length(h1),mu(2)),'.g','markersize',5)
t = legend('X_1','X_1 Poisson','X_2','X_2 Poisson');
set(t,'box','off')
axis square
xlabel('X'), ylabel('P(X)')
title(sprintf('Step 4: Correlated Poisson, Positive Correlation\nMarginal Distributions with Poisson\n distribution with correct mean'))

% Finally, we convince ourselves that the two variables are indeed
% positively correlated by looking at the 2D joint histogram. We can see
% that although the marginals are perfect Poissonian, most samples are
% concentrated along the main diagonal, indicating the positive
% correlation.

subplot(2,2,3)
hh = joints2D{2};   % joint histogram
imagesc(hh)
axis square
title('')
xlabel('X_1'), ylabel('X_2')

disp('To proceed to negatively correlated Poisson variables, hit any key...')
pause

%% Step 5: Generating correlated Poisson variables via the method
% in section 3.3 for negative correlations
%
% Again, we generate samples from a correlated Poisson distribution with
% specified mean and covariance matrix. This time the covariance is
% negative though.

fprintf('\nStep 5: Generating negatively correlated Poisson variables (section 3.3)\n\n')
fprintf('Computing...\n')

mu = [7,9]';      % set the mean
C = [7 -3;-3 9];    % set an admissable covariance matrix

[S,gamma,Lambda,joints2D] = DGPoisson(mu,C,1e5);   % generate sample via discretized Gaussian

% Next, we verify that the samples have desired mean and covariance
% structure. We see that this is indeed the case.

muHat = mean(S,2);
CHat = cov(S');

fprintf('Target mean X1:  %.2f     Estimated mean X1:  %.3f\n',mu(1),muHat(1))
fprintf('Target mean X2:  %.2f     Estimated mean X2:  %.3f\n',mu(2),muHat(2))
fprintf('Target cov X1X2: %.2f     Estimated cov X1X2: %.3f\n',C(1,2),CHat(1,2))

disp('To plot the marginals and joint histogram, hit any key...')
pause

figure(f)

% Again we compare the marginal distributions to the Poisson distribution
% with the same mean and we find them matching very well.

subplot(2,2,2)
h1 = joints2D{1,1};   % marginal of X_1
h2 = joints2D{2,2};   % marginal of X_2
plot(0:length(h1)-1,h1,'-r','markersize',5), hold on
plot(0:length(h1),poisspdf(0:length(h1),mu(1)),'.r','markersize',5)
plot(0:length(h2)-1,h2,'-g','markersize',5)
plot(0:length(h1),poisspdf(0:length(h1),mu(2)),'.g','markersize',5)
t = legend('X_1','X_1 Poisson','X_2','X_2 Poisson');
set(t,'box','off')
axis square
xlabel('X'), ylabel('P(X)')
title(sprintf('Step 5: Correlated Poisson, Negative Correlation\nMarginal Distributions with Poisson\n distribution with correct mean'))

% Nevertheless, the samples have a different structure as we can see from
% the 2D joint histogram. This time, the two variables are negatively
% correlated and when X_1 tends to take large values, X_2 tends to take low
% ones.

subplot(2,2,4)
hh = joints2D{2};   % joint histogram
imagesc(hh)
axis square
title('')
xlabel('X_1'), ylabel('X_2')

author: Sable

repo: mcbench-benchmarks


text:
%Program for SVD-based Image Quality Measure

%Program Description
% This program outputs a graphical & numerical image quality measure
% based on Singular Value Decomposition.
% For details on the implementation, please refer
% Aleksandr Shnayderman, Alexander Gusev, and Ahmet M. Eskicioglu,
% "An SVD-Based Grayscale Image Quality Measure for Local and Global Assessment",
% IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 15, NO. 2, FEBRUARY 2006.
%
%Parameters
% refImg        -   Input Reference Gray Image
% distImg       -   Input Distorted Gray Image
% blkSize       -   Window size for block processing
% graMeasure    -   Graphical Image quality measure
% scaMeasure    -   Numerical Image quality measure
%
%Author : Athi Narayanan S
%Student, M.E, EST,
%K.S.R College of Engineering
%Erode, Tamil Nadu, India.
%s_athi1983@yahoo.co.in
%http://sites.google.com/site/athisnarayanan/

function [graMeasure, scaMeasure] = SVDQualityMeasure(refImg, distImg, blkSize)

k = size(refImg, 1);

blkx = blkSize;
blky = blkSize;

blockwise1 = MatDec(refImg,blkx);
blockwise2 = MatDec(distImg,blkx);
[blkx blky imgx imgy] = size(blockwise1);
graMeasure = zeros(imgx,imgy);
blockwise1 = double(blockwise1);
blockwise2 = double(blockwise2);

for i=1:imgx
    for j=1:imgy
        temp_in_image = blockwise1(:,:,i,j);temp_in_image=temp_in_image(:);
        original_img = reshape(temp_in_image,blkx,blky);
        temp_dist_image = blockwise2(:,:,i,j);temp_dist_image=temp_dist_image(:);
        distorted_img = reshape(temp_dist_image,blkx,blky);
        graMeasure(i,j) = sqrt(sum((svd(original_img)-svd(distorted_img)).^2));
    end
end

graMeasure = round((graMeasure/max(max(graMeasure)))*255);
graMeasure = uint8(graMeasure);

scaMeasure = sum(sum(abs(graMeasure-median(median(graMeasure)))))/((k/blkx).^2);