text:
% nsamples  calculate the number of samples yet needed
%
% N = nsamples(no_i,ptNum,s,conf)
% no_i ... current number of inliers
% ptNum ... total number of points
% s ... sample size
% conf ... confidence value
%
% $Id: nsamples.m,v 1.1 2005/05/23 16:15:59 svoboda Exp $

function N = nsamples(no_i,ptNum,s,conf)

outl = 1-no_i/ptNum;
N	 = log(1-conf) / log(1-(1-outl)^s+eps);

return;

author: strawlab

repo: MultiCamSelfCal


text:
function output = F_sparse_affine_transform(input, transform, bias, singlePrecision)

[D,M,N] = size(input);
if N>1
    input = reshape(input, D,M*N);
end

visible_nonzero_idx = find(sum(abs(input),2)>0);
visible_nonzero = full(input(visible_nonzero_idx,:));
if singlePrecision==1
    visible_nonzero = single(visible_nonzero);
end
output = bsxfun(@plus, transform(:,visible_nonzero_idx) * visible_nonzero, bias);    % only compute the nonzero elements

if N>1
    output = reshape(output, size(output,1), M,N);
end

end

author: singaxiong

repo: SignalGraph


text:
classdef (TestTags = {'Unit', 'SPGR', 'qMT'}) BlochSol_Test < matlab.unittest.TestCase
%% BLOCHSOL_TEST Test class for BlochSol.m
%
%   --tests--
%   test_no_pulse_mag_remains_unchanged
%       - Assert that a magnetization vector with only longitudinal values
%         remains constant for a constant pulse of 0 degree FA.
%
%   test_no_longi_mag_evol_for_transverse_initial_mag_inf_t2_t1
%       - Assert that a magnetization vector with only transverse values
%         does not accrue longitudinal magnetization for the case of
%         infinite T1 and T2 (0 deg FA).
%
%   test_magn_evols_to_equilibrium_vals
%       - Assert that an empty magnetization vector evolves to the
%         equilibrium longitudinal values (M0f and M0r) after a time
%         much greater than T1f (10x) and 0 deg FA.
%
%   test_mag_goes_to_zero_for_time_much_greater_than_t2f
%       - Assert that a magnetization vector with only transverse values
%         decays to 0 after a time much greater than T2f (x10)
%
%   test_steady_state_smaller_for_small_delta_than_big
%       - Assert that steady-state longitudinal (free) magnization is
%         smaller for the case of a small offset frequency than large
%         offset freq.
%
%   test_steady_state_larger_for_small_FA_than_big
%       - Assert that steady-state longitudinal (free) magnization is
%         larger for the case of a small flip angle (per millisec) than
%         large flip angle.
%

    properties
        Param = struct('M0f', 1,           ...
                       'M0r', 0.15,        ...
                       'R1f', 1.1,         ...
                       'R1r', 1,           ...
                       'R2f', 1/0.03,      ...
                       'R2r', 1/(12*10^-6),...
                       'kf',  4.0,         ...
                       'kr',  4.0/0.15);
    end

    methods (TestClassSetup)
    end

    methods (TestClassTeardown)
    end

    methods (Test)
        function test_no_pulse_mag_remains_unchanged(testCase)
            %% Prep
            %
            delta = 200; % On-resonance
            testCase.Param.G = computeG(delta,1/testCase.Param.R2f, 'SuperLorentzian');

            pulseDur = 10; % seconds
            pulseShape = 'hard';

            flipAngle = 0; % No pulse
            Pulse = GetPulse(flipAngle, delta, pulseDur, pulseShape);

            M0 = [0 0 1 0.15]'; % Initial magnetization, Xf Yf Zf Zr

            timeRange = linspace(0,10,1000); % in seconds

            for ii =1:length(timeRange)
               M(ii, :) = BlochSol(timeRange(ii), M0, testCase.Param, Pulse);
            end

            %                    Actual ,  Expected
            testCase.assertEqual(M(end,:), M0', 'AbsTol', 0.001);
        end

        function test_no_longi_mag_evol_for_transverse_initial_mag_inf_t2_t1(testCase)
            %% Prep
            %
            Param = testCase.Param;

            delta = 100; % Go to the rotating ref frame for this case
            Param.G = computeG(delta,1/Param.R2f, 'SuperLorentzian');

            pulseDur = 0.001; % seconds, near-instantaneous
            pulseShape = 'hard';

            flipAngle = 0;
            Pulse = GetPulse(flipAngle, delta, pulseDur, pulseShape);

            %***
            Param.R2f = 0.0001; % Assume very large T2f
            Param.R2r = 0.0001; % Assume very large T2r
            Param.R1f = 0.0001; % Assume very large T1f
            Param.R1r = 0.0001; % Assume very large T1r
            %***

            M0 = [1 0 0 0]'; % Initial magnetization, Xf Yf Zf Zr

            timeRange = linspace(0, 1, 1000); % in seconds

            for ii =1:length(timeRange)
               M(ii, :) = BlochSol(timeRange(ii), M0, Param, Pulse);
            end

            %                    Actual                     , Expected
            testCase.assertEqual([M(end,3) M(end,4)], [M0(3) M0(4)], 'AbsTol', 0.001);
        end

        function test_magn_evols_to_equilibrium_vals(testCase)
            %% Prep
            %
            Param = testCase.Param;

            delta = 100; % Go to the rotating ref frame for this case
            Param.G = computeG(delta,1/Param.R2f, 'SuperLorentzian');

            pulseDur = 0.001; % seconds, near-instantaneous
            pulseShape = 'hard';

            flipAngle = 0;
            Pulse = GetPulse(flipAngle, delta, pulseDur, pulseShape);

            M0 = [0 0 0 0]'; % Initial magnetization, Xf Yf Zf Zr

            timeRange = linspace(0, (1/Param.R1f*10), 1000); % in seconds, 10x T1f

            for ii =1:length(timeRange)
               M(ii, :) = BlochSol(timeRange(ii), M0, Param, Pulse);
            end

            %                    Actual ,  Expected
            testCase.assertEqual(M(end,:), [0 0 Param.M0f Param.M0r], 'AbsTol', 0.001);
        end

        function test_mag_goes_to_zero_for_time_much_greater_than_t2f(testCase)
            %% Prep
            %
            Param = testCase.Param;

            delta = 100; % Go to the rotating ref frame for this case
            Param.G = computeG(delta,1/Param.R2f, 'SuperLorentzian');

            pulseDur = 0.001; % seconds, near-instantaneous
            pulseShape = 'hard';

            flipAngle = 0;
            Pulse = GetPulse(flipAngle, delta, pulseDur, pulseShape);

            M0 = [1 0 0 0]'; % Initial magnetization, Xf Yf Zf Zr

            timeRange = linspace(0, (1/Param.R2f)*10, 1000); % in seconds, 10x T2f

            for ii =1:length(timeRange)
               M(ii, :) = BlochSol(timeRange(ii), M0, Param, Pulse);
            end

            %                    Actual ,  Expected
            testCase.assertEqual(M(end,1:2), [0 0], 'AbsTol', 0.001);
        end

        function test_steady_state_smaller_for_small_delta_than_big(testCase)
            %% Prep
            %
            Param1 = testCase.Param;
            Param2 = testCase.Param;

            delta1 = 100;
            delta2 = 10*delta1;

            Param1.G = computeG(delta1,1/Param1.R2f, 'SuperLorentzian');
            Param2.G = computeG(delta2,1/Param2.R2f, 'SuperLorentzian');

            pulseDur = 10; % seconds
            pulseShape = 'hard';

            flipAngle = 40; % FA per millisec
            Pulse1 = GetPulse((flipAngle*1000)*pulseDur, delta1, pulseDur, pulseShape);
            Pulse2 = GetPulse((flipAngle*1000)*pulseDur, delta2, pulseDur, pulseShape);

            M0 = [0 0 testCase.Param.M0f testCase.Param.M0r]'; % Initial magnetization, Xf Yf Zf Zr

            timeRange = linspace(0, pulseDur, 1000); % in seconds, 10x T2f

            for ii =1:length(timeRange)
               M1(ii, :) = BlochSol(timeRange(ii), M0, Param1, Pulse1);
               M2(ii, :) = BlochSol(timeRange(ii), M0, Param2, Pulse2);
            end

            %                          Small delta, Big delta
            testCase.assertLessThan(M1(end,3)  , M2(end,3));
        end

        function test_steady_state_larger_for_small_FA_than_big(testCase)
            %% Prep
            %
            Param = testCase.Param;

            delta = 1000;

            Param.G = computeG(delta,1/Param.R2f, 'SuperLorentzian');

            pulseDur = 10; % seconds
            pulseShape = 'hard';

            flipAngle1 = 40; % FA per millisec
            flipAngle2 = 2*flipAngle1 ; % FA per millisec

            Pulse1 = GetPulse((flipAngle1*1000)*pulseDur, delta, pulseDur, pulseShape);
            Pulse2 = GetPulse((flipAngle2*1000)*pulseDur, delta, pulseDur, pulseShape);

            M0 = [0 0 testCase.Param.M0f testCase.Param.M0r]'; % Initial magnetization, Xf Yf Zf Zr

            timeRange = linspace(0, pulseDur, 1000); % in seconds, 10x T2f

            for ii =1:length(timeRange)
               M1(ii, :) = BlochSol(timeRange(ii), M0, Param, Pulse1);
               M2(ii, :) = BlochSol(timeRange(ii), M0, Param, Pulse2);
            end

            %                          Small FA , Big FA
            testCase.assertGreaterThan(M1(end,3), M2(end,3));
        end
    end

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [G] = spm_lx_phase (P,M)
% Observation function for phase-coupled oscillators
% FORMAT [G] = spm_lx_phase (P,M)
%
% G     Observations y = Gx

Nr=length(P.L);
G=eye(Nr);

author: spm


text:
% cmbglobl.m : cmb globals : physical constants that are used through all files
% are given a value here
%
% D Vangheluwe 28 jan 2005, see cmbaccur.m anc cmbacc1.m as the main files

global  GL_cmb_c GL_cmb_h0  GL_cmb_t0  GL_cmb_T0  GL_cmb_rv  GL_cmb_fv  GL_cmb_kg1  GL_cmb_ka1 ...
  GL_cmb_yp  GL_cmb_ncr  GL_cmb_cr  GL_cmb_pcm  GL_cmb_dha;

% velocity of light in m/s
GL_cmb_c = 2.998e8;
% the Hubble constant at present in (h Mpc^-1), see my notes p71
GL_cmb_h0 = 1e5/GL_cmb_c;
% the present temperature of the CMB in degr Kelvin and in eV
GL_cmb_t0 = 2.725;
GL_cmb_T0 = GL_cmb_t0/11605;
% ratio of radiation density/critical density, see Dodelson (2.87)
%fv = 0.405;
GL_cmb_rv = (21/8) * (4/11)^(4/3);
GL_cmb_fv = GL_cmb_rv/(1 + GL_cmb_rv);
GL_cmb_kg1 = 2.47e-5;
GL_cmb_ka1 = GL_cmb_kg1/(1 - GL_cmb_fv);
%ka1 = 4.15e-5;
% the primordial helium mass fraction
GL_cmb_yp = 0.24;
% critical density in m^-3
GL_cmb_ncr = 11.23;
%ncr = 10.8;
% Compton cross section in m^2, see Dodelson p 72
GL_cmb_cr = 0.665e-28;
% conversion of Parsec to meters
GL_cmb_pcm = 3.0856e22;
% amplitude of the primordial density perturbation at the horizon, see Dodelson (8.76)
GL_cmb_dha = 4.47e-5;

author: Sable

repo: mcbench-benchmarks


text:
function [forecast_record]=maforecast(data_endo_a,data_exo_a,data_exo_p,It,Bu,beta_gibbs,sigma_gibbs,delta_gibbs,Fperiods,n,m,p,k1,k3)

% function [forecast_record]=maforecast(data_endo_a,data_exo_a,data_exo_p,It,Bu,beta_gibbs,sigma_gibbs,delta_gibbs,Fperiods,n,m,p,k1,k3)
% computes draws from the posterior predictive distribution, that is, from the posterior distribution of forecasts
% inputs:  - matrix 'data_endo_a': the matrix storing the pre-forecast endogenous data
%          - matrix 'data_exo_a': the matrix storing the pre-forecast exogenous data
%          - matrix 'data_exo_p': the matrix storing the predicted exogenous data
%          - integer 'It': the total number of iterations run by the Gibbs sampler
%          - integer 'Bu': the number of initial iterations discared as burn-in sample
%          - matrix 'beta_gibbs': the matrix recording the post-burn draws of beta
%          - matrix 'sigma_gibbs': the matrix recording the post-burn draws of sigma
%          - matrix 'delta_gibbs': the matrix recording the post-burn draws of delta
%          - integer 'Fperiods': the number of periods for which forecasts have to be produced
%          - integer 'n': the number of endogenous variables in the model
%          - integer 'm': the number of exogenous variables in the model
%          - integer 'p': the number of lags in the model
%          - integer 'k1': the number of coefficients related to the endogenous variables for each equation in the model
%          - integer 'k3': the number of coefficients related to the exogenous variables for each equation, in the reformulated model (3.5.5)
% outputs: - cell 'forecast_record': the cell array containing records of simulated  forecasts

% this function implements algorithm 2.1.1, adapted for the mean-adjusted BVAR model

% create first the cell storing the results
forecast_record=cell(n,1);

% generate the matrix of predicted exogenous variables
% augment the matrices of exogenous with a column of ones to account for the exogenous
data_exo_a=[ones(size(data_endo_a,1),1) data_exo_a];
data_exo_p=[ones(Fperiods,1) data_exo_p];

% then start simulations
% repeat the process a number of times equal to the number of simulations retained from Gibbs sampling
for ii=1:It-Bu

% compute the matrix temp1
temp1=data_endo_a(end-p+1:end,:);

% compute the matrix temp2
temp2=[data_exo_a(end-p+1:end,:)];

% step 2: draw beta and sigma
% draw beta from its posterior distribution
beta=beta_gibbs(:,ii);
% reshape
B=reshape(beta,k1,n);

% draw delta from its posterior distribution
delta=delta_gibbs(:,ii);
% reshape
Delta=reshape(delta,k3,n);

% draw sigma from its posterior distribution
sigma=sigma_gibbs(:,ii);
% reshape sigma to recover its original square form
sigma=reshape(sigma,n,n);

   % step 4: generate forecasts recursively
   % repeat the process for periods T+1 to T+h
   for jj=1:Fperiods

   % concatenate the predicted exogenous to the top of temp2 (the actual exogenous)
   temp2=[temp2;data_exo_p(jj,:)];

   % use the function lagx on temp1 to obtain the matrix Y; retain only the last row
   X=bear.lagx(temp1,p-1);
   X=X(end,:);

   % use the function lagx on temp2 to obtain the matrix Z; retain only the last row
   Z=bear.lagx(temp2,p);
   Z(:,m+1:end)=-Z(:,m+1:end);
   Z=Z(end,:);

   % draw the residuals from N(0,sigma)
   res=bear.trns(chol(bear.nspd(sigma),'Lower')*randn(n,1));

   % obtain predicted value for T+jj by using (3.5.9)
   yp=X*B+Z*Delta+res;

   % concatenate the transpose of yp to the top of temp1
   temp1=[temp1;yp];

   % repeat until values are obtained for T+h
   end

% step 5: record the results from current iteration in the cell forecast_record
   % loop over variables
   for kk=1:n
   % consider column kk of matrix temp1 and select the last h rows: these are the predicted values for the period T+1 to T+h, for variable kk
   temp3=temp1(end-Fperiods+1:end,kk);
   % record these values in the corresponding matrix of forecast_record
   forecast_record{kk,1}(ii,:)=temp3';
   end

% then go for next iteration
end

author: european-central-bank

repo: BEAR-toolbox


text:
function plot_average_frequency(w)
    figure
    sp=spectralobject(1024,924,100,[]);
    % Plot a spectrogram with superposed frequency metrics
    [Tcell,meanF,peakF]=spectrogram(w,'spectralobject',sp,'plot_metrics',1);
    % Plot the frequency metrics on their own
    for c=1:numel(w)
        subplot(numel(w),1,c)
        plot(Tcell{c},smooth(peakF{c}),'g')
        hold on
        plot(Tcell{c},meanF{c},'k')
        datetick('x')
        ylabel('Hz')
        sta= get(w(c),'station');
        chan= get(w(c),'channel');
        th=text(0.1,0.9, sprintf('%s %s.%s',datestr(Tcell{c}(1),30),sta,chan),'Units','normalized')
    end

end

author: geoscience-community-codes

repo: GISMO


text:
function [ n_data, h, a, t ] = owen_values ( n_data )

%*****************************************************************************80
%
%% OWEN_VALUES returns some values of Owen's T function.
%
%  Discussion:
%
%    Owen's T function is useful for computation of the bivariate normal
%    distribution and the distribution of a skewed normal distribution.
%
%    Although it was originally formulated in terms of the bivariate
%    normal function, the function can be defined more directly as
%
%      T(H,A) = 1 / ( 2 * pi ) *
%        Integral ( 0 <= X <= A ) e^(H^2*(1+X^2)/2) / (1+X^2) dX
%
%    In Mathematica, the function can be evaluated by:
%
%      fx = 1/(2*Pi) * Integrate [ E^(-h^2*(1+x^2)/2)/(1+x^2), {x,0,a} ]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 May 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Mike Patefield, David Tandy,
%    Fast and Accurate Calculation of Owen's T Function,
%    Journal of Statistical Software,
%    Volume 5, Number 5, 2000, pages 1-25.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real H, a parameter.
%
%    Output, real A, the upper limit of the integral.
%
%    Output, real T, the value of the function.
%
  n_max = 28;

  a_vec = [ ...
    0.2500000000000000E+00, ...
    0.4375000000000000E+00, ...
    0.9687500000000000E+00, ...
    0.0625000000000000E+00, ...
    0.5000000000000000E+00, ...
    0.9999975000000000E+00, ...
    0.5000000000000000E+00, ...
    0.1000000000000000E+01, ...
    0.2000000000000000E+01, ...
    0.3000000000000000E+01, ...
    0.5000000000000000E+00, ...
    0.1000000000000000E+01, ...
    0.2000000000000000E+01, ...
    0.3000000000000000E+01, ...
    0.5000000000000000E+00, ...
    0.1000000000000000E+01, ...
    0.2000000000000000E+01, ...
    0.3000000000000000E+01, ...
    0.5000000000000000E+00, ...
    0.1000000000000000E+01, ...
    0.2000000000000000E+01, ...
    0.3000000000000000E+01, ...
    0.5000000000000000E+00, ...
    0.1000000000000000E+01, ...
    0.2000000000000000E+01, ...
    0.3000000000000000E+01, ...
    0.1000000000000000E+02, ...
    0.1000000000000000E+03 ];

  h_vec = [ ...
    0.0625000000000000E+00, ...
    6.5000000000000000E+00, ...
    7.0000000000000000E+00, ...
    4.7812500000000000E+00, ...
    2.0000000000000000E+00, ...
    1.0000000000000000E+00, ...
    0.1000000000000000E+01, ...
    0.1000000000000000E+01, ...
    0.1000000000000000E+01, ...
    0.1000000000000000E+01, ...
    0.5000000000000000E+00, ...
    0.5000000000000000E+00, ...
    0.5000000000000000E+00, ...
    0.5000000000000000E+00, ...
    0.2500000000000000E+00, ...
    0.2500000000000000E+00, ...
    0.2500000000000000E+00, ...
    0.2500000000000000E+00, ...
    0.1250000000000000E+00, ...
    0.1250000000000000E+00, ...
    0.1250000000000000E+00, ...
    0.1250000000000000E+00, ...
    0.7812500000000000E-02, ...
    0.7812500000000000E-02, ...
    0.7812500000000000E-02, ...
    0.7812500000000000E-02, ...
    0.7812500000000000E-02, ...
    0.7812500000000000E-02 ];

  t_vec = [ ...
    3.8911930234701366E-02, ...
    2.0005773048508315E-11, ...
    6.3990627193898685E-13, ...
    1.0632974804687463E-07, ...
    8.6250779855215071E-03, ...
    6.6741808978228592E-02, ...
    0.4306469112078537E-01, ...
    0.6674188216570097E-01, ...
    0.7846818699308410E-01, ...
    0.7929950474887259E-01, ...
    0.6448860284750376E-01, ...
    0.1066710629614485E+00, ...
    0.1415806036539784E+00, ...
    0.1510840430760184E+00, ...
    0.7134663382271778E-01, ...
    0.1201285306350883E+00, ...
    0.1666128410939293E+00, ...
    0.1847501847929859E+00, ...
    0.7317273327500385E-01, ...
    0.1237630544953746E+00, ...
    0.1737438887583106E+00, ...
    0.1951190307092811E+00, ...
    0.7378938035365546E-01, ...
    0.1249951430754052E+00, ...
    0.1761984774738108E+00, ...
    0.1987772386442824E+00, ...
    0.2340886964802671E+00, ...
    0.2479460829231492E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    h = 0.0;
    a = 0.0;
    t = 0.0;
  else
    h = h_vec(n_data);
    a = a_vec(n_data);
    t = t_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function XML=mat2xml(MAT,VARNAME)

% MAT2XML converts structured variable MAT into XML string
%
%Syntax XML=mat2xml(MAT,VARNAME)
%
%Description
%  MAT : structured varable
%  VARNAME : variable name (string)
%  XML : xml version of structured variable (string)
%
% See Also: XML2MAT
%
% Jonas Almeida, almeidaj@musc.edu, 20 Aug 2002, XML4MAT Tbox

if nargin<2;VARNAME='ans';end % if not provided make it a matlab answer variable
w=whos('MAT');w.name=VARNAME;

XML=['<',w.name,' class="',w.class,'" size="',num2str(w.size),'">'];
if strcmp(w.class,'char')
    XML=[XML,spcharin(MAT(:)')];
elseif strcmp(w.class,'struct')
    names=fieldnames(MAT);
    %struct_fields=[' fields="',names{1}];for j=2:length(names);struct_fields=[struct_fields,' ',names{j}];end;struct_fields=[struct_fields,'">'];XML=[XML(1:(end-1)),struct_fields];
    for i=1:prod(w.size)
        for j=1:length(names)
            XML=[XML,mat2xml(eval(['MAT(i).',names{j}]),names{j})];
        end
    end
elseif strcmp(w.class,'cell')
    for i=1:prod(w.size)
        XML=[XML,mat2xml(MAT{i},'cell')];
    end
else %if strcmp(w.class,'double')|strcmp(w.class,'single')
    XML=[XML,num2str(MAT(:)')];
end
XML=[XML,'</',w.name,'>'];
XML=regexprep(XML,'\s{2,}',' ');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%-----------------------------------------------------------------------
% Job saved on 22-Jul-2014 14:45:14 by cfg_util (rev $Rev$)
% spm SPM - SPM12b (5672)
% cfg_basicio BasicIO - Unknown
%-----------------------------------------------------------------------
matlabbatch{1}.spm.spatial.normalise.write.subj.def = {'/Users/kasperla/Documents/code/matlab/fmri_svn/development/analysis/examples/model_based_fmri_3T/zFatTmp/y_raw.nii'};
matlabbatch{1}.spm.spatial.normalise.write.subj.resample = {'/Users/kasperla/Documents/code/matlab/fmri_svn/development/analysis/examples/model_based_fmri_3T/zFatTmp/mraw.nii,1'};
matlabbatch{1}.spm.spatial.normalise.write.woptions.bb = [-78 -112 -70
                                                          78 76 85];
matlabbatch{1}.spm.spatial.normalise.write.woptions.vox = [2 2 2];
matlabbatch{1}.spm.spatial.normalise.write.woptions.interp = 7;

author: translationalneuromodeling

repo: tapas


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% TLD is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with TLD.  If not, see <http://www.gnu.org/licenses/>.

% Source for TLD demo (CVPR'10)
% Usage:
%  1. Define input
%     a) plug in your camera
%     b) fill '/_input' with image sequence, define 'init.txt' (initial bounding box position)
%  2. run demo.m
%  3. Output of the TLD is written to 'log.txt', format [frame id; left col; top row; right col; bottom row]
%
%
% Control keys:
%    n ... shows negative examples in online model (default on)
%    p ... shows positive examples in online model (default on)
%    i ... initialization of different target
%    c ... show confidence score (default on)
%    o ... show output as circle/dot/no output (default circle)
%    d ... show detections (default on)
%    t ... show target in top left corner (default off)
%    r ... replace target with first patch (default off)
%    # ... draw trajectory of target (default off)
%    <space> ... svae current image
%    1 ... mode without learning (fast, does not improve)
%    2 ... mode with learning (slower, but can learn new appearances and discriminate)
%    q ... quit application

addpath(genpath('.')); init_workspace;

% Input
opt.source          = struct('camera',0,'input','_input/','bb0',[]);
opt.debug           = 0;
opt.name            = '_snapshots/'; mkdir(opt.name); % saves s snapshot every 100th frames

% Parameters
opt.model           = struct('ncc_thesame',0.95,'valid',0.5,'patchsize',[10 10],'min_win',20,'num_trees',10,'num_features',13,'thr_fern',0.5,'thr_nn',0.65,'thr_nn_valid',0.7,'fliplr',0);
opt.p_par_init      = struct('num_closest',10,'num_warps',20,'noise',5,'angle',20,'shift',0.02,'scale',0.02);
opt.p_par_update    = struct('num_closest',10,'num_warps',10,'noise',5,'angle',10,'shift',0.02,'scale',0.02);
opt.n_par           = struct('overlap',0.2,'num_patches',100,'num_synthetic',0,'noise',1);
opt.tracker         = struct('occlusion',10,'grid',10,'big_fb',10);
opt.control         = struct('maxbbox',.9,'update_detector',1,'drop_img',1,'repeat',1,'rescale',1);
opt.plot            = struct('pex',1,'nex',1,'save',0,'dt',0,'confidence',0,'target',0,'replace',0,'drawoutput',3,'draw',0,'pts',0,'help', 0,'patch_rescale',1);

%  profile on;
global tld;
tldDemo(opt);
%   profile off;
%   profile viewer;

repo: OpenTLD


text:
function E=fit2(x0,x,y)
E=sum(abs( x0(1)*x+x0(2)-y ));

author: dynamicslab


text:
% std_moveoutlier()  - Commandline function, to reassign specified outlier component(s)
%                            from a cluster to its outlier cluster.
% Usage:
%                   >> STUDY = std_moveoutlier(STUDY, ALLEEG, from_cluster, comps);
% Inputs:
%   STUDY         - EEGLAB STUDY set comprising some or all of the EEG datasets in ALLEEG.
%   ALLEEG        - global EEGLAB vector of EEG structures for the dataset(s) included in the STUDY.
%                       ALLEEG for a STUDY set is typically created using load_ALLEEG().
%   from_cluster - cluster number, the cluster outlier components are moved from.
%   comps          - [numeric vector] component indices in the from_cluster to move.
%
% Outputs:
%   STUDY    - the input STUDY set structure modified with the components reassignment.
%
%   Example:
%                         >> from_cluster = 10; comps = [2 7];
%                         >> STUDY = std_movecomp(STUDY,ALLEEG, from_cluster, to_cluster, comps);
%                    Components 2 and 7 of cluster 10 are moved to the its outlier cluster ('Outliers Cls 10').
%
%  See also  pop_clustedit
%
% Authors:  Hilit Serby, Arnaud Delorme, Scott Makeig, SCCN, INC, UCSD, June, 2005

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

function STUDY = std_moveoutlier(STUDY, ALLEEG, old_clus, comps)

% Cannot move components if the cluster is either a 'Notclust' or
% 'Outliers' cluster
if strncmpi('Notclust',STUDY.cluster(old_clus).name,8) | strncmpi('Outliers',STUDY.cluster(old_clus).name,8)
    warndlg2('std_moveoutlier: cannot move components from Notclust or Outliers cluster');
    return;
end

% Cannot move components if clusters have children clusters
if ~isempty(STUDY.cluster(old_clus).child)
    warndlg2('Cannot move components if cluster has children clusters!' , 'Aborting remove outliers');
    return;
end

outlier_clust = std_findoutlierclust(STUDY,old_clus); %find the outlier cluster for this cluster
if outlier_clust == 0 %no such cluster exist
    STUDY = std_createclust(STUDY, ALLEEG, ['Outliers ' STUDY.cluster(old_clus).name]); %create an outlier cluster
    outlier_clust = length(STUDY.cluster);
end

%move the compnents to the outliers cluster
STUDY = std_movecomp(STUDY, ALLEEG, old_clus, outlier_clust, comps);

author: PatternRecognition

repo: OpenBMI


text:
function [At,b,Z] = getequation(symexpr,vartable,decvartable,varmat)
%function [At,b,Z] = getequation(symexpr,vartable,decvartable,decvartablename)
%
% GETEQUATION --- Convert a symbolic expression to At, b, and Z
%         used in an SOS program.

% This file is part of SOSTOOLS - Sum of Squares Toolbox ver 3.00.
%
% Copyright (C)2002, 2004, 2013  A. Papachristodoulou (1), J. Anderson (1),
%                                G. Valmorbida (1), S. Prajna (2),
%                                P. Seiler (3), P. A. Parrilo (4)
% (1) Department of Engineering Science, University of Oxford, Oxford, U.K.
% (2) Control and Dynamical Systems - California Institute of Technology,
%     Pasadena, CA 91125, USA.
% (3) Aerospace and Engineering Mechanics Department, University of
%     Minnesota, Minneapolis, MN 55455-0153, USA.
% (4) Laboratory for Information and Decision Systems, M.I.T.,
%     Massachusetts, MA 02139-4307
%
% Send bug reports and feedback to: sostools@cds.caltech.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%

% 03/01/02 - SP
% 03/10/02 - SP -- Use Maple
% 07/10/13 - JA&GV -- Use the matlab's symbolic engine
% 09/11/13 - PJS -- Addition for multipoly objects

if isa(symexpr,'polynomial')
    % PJS: Handle Polynomial Objects including the Matrix Case
    decvarnum = numel(decvartable);
    vartable = [vartable; varmat];
    cvartable = char(vartable);
    dimp = size(symexpr,2);

    % Collect symexpr = g(c)*h(x) where x are the vars in vartable,
    % c are the decision variables, and h(x) is a vector of monoms.
    % in the actual polynomial vars.  Use this to find unique
    % monomials in the polynomial variables.
    [g0,g,h] = collect(symexpr(:),setdiff(symexpr.varname,cvartable));
    g = g(:);
    if ~isequal(g0,0)
        g = [g0;g];
        h = [1; h];
    end
    [nmon,nvar] = size(h.degmat);

    % Reorder the monomials matrix with variables in the order
    % listed in cvartable and sorted monomials
    Z = zeros(nmon,nvar);
    [~,idx]=ismember(h.varname,cvartable);
    Z(:,idx) = full(h.degmat);
    Z = sortNoRepeat([],Z);

    % Initialize coefficient variables
    [nmon,nvar] = size(Z);
    coeffnts = sparse(nmon*dimp,dimp);
    if decvarnum>0
        coeffnts_decvar = polynomial(sparse(nmon*dimp,dimp));
    end

    % Define coefficients
    for i = 1:dimp
        for j = i:dimp
            exprij = symexpr(i,j);
            [g0,g,h] = collect(exprij,setdiff(exprij.varname,cvartable));
            g = g(:);
            if ~isequal(g0,0)
                g = [g0;g];
                h = [1; h];
            end
            coefmatr = double(subs(g,decvartable,zeros(decvarnum,1)));
            monmatr = zeros(length(h),nvar);
            [~,idx]=ismember(h.varname,cvartable);
            monmatr(:,idx) = h.coef'*h.degmat;

            for k = 1:length(h);
                s_ijk = coefmatr(k,1);
                s_ijk_decvar = g(k);
                mon_k= monmatr(k,:);

                [val,ind_k] = max(sum((Z == kron(ones(nmon,1),mon_k)),2));
                coeffnts((ind_k-1)*dimp+i,j) = s_ijk;
                coeffnts((ind_k-1)*dimp+j,i) = s_ijk;

                coeffnts_decvar((ind_k-1)*dimp+i,j) = s_ijk_decvar;
                coeffnts_decvar((ind_k-1)*dimp+j,i) = s_ijk_decvar;
            end
        end
    end

    % Constructing At and b matrices
    At = sparse(decvarnum,size(Z,1)*dimp^2);
    if decvarnum>0
        for i=1:nmon
            Mivec = reshape(coeffnts_decvar((i-1)*dimp+1:i*dimp,:),dimp^2,1);
            At(:,(i-1)*dimp^2+1:i*dimp^2) = sparse(-double(jacobian(Mivec,decvartable))');
        end
    end
    b = sparse(coeffnts);

else
    % PJS: Original Code to Handle Symbolic Objects

    if strcmp(decvartable,'[]');
        decvarnum = 0;
    else
        decvarnum = length(find(decvartable == ','))+1;
    end;
    expr = evalin(symengine,symexpr);
    %vartable = evalin(symengine,vartable);
    if  length(varmat)==2
        vartable = evalin(symengine,vartable);
    else
        vartable = evalin(symengine,[vartable(1:end-1),',',varmat(2:end)]);  %JA&GV 07/06/13
    end

    decvartable = evalin(symengine,decvartable);
    charvartable = converttochar([vartable]);

    FPexpr = feval(symengine,'expand',expr);

    %JA&GV 07/13 commented the below to implement the matrix case
    %FPexpr = feval(symengine,'collect',FPexpr,charvartable);

    dimp = size(FPexpr,2);%JA&GV jul/04/2013 sets the dimension of the matrix of polynomials

    for i = 1:size(FPexpr,1)
        for j = 1:dimp
            FPexpr(i,j) = feval(symengine,'collect',FPexpr(i,j),charvartable);
        end
    end

    if isempty(decvartable)==0%JA&GV 07/13 the code below  addresses the matrix case

        Zfull = [];
        for i = 1:dimp
            for j = i:dimp
                coefmon = feval(symengine,'poly2list',FPexpr(i,j),charvartable);
                coefmonmatr = subs(coefmon.',decvartable,zeros(1,length(decvartable)));
                Z = double(coefmonmatr(:,2:end));
                Zfull = sortNoRepeat(Zfull,Z);
            end
        end
        Z = Zfull;
        [nmon,nvar] = size(Z);
        coeffnts = sparse(nmon*dimp,dimp);
        coeffnts_decvar = sym(sparse(nmon*dimp,dimp));

        for i = 1:dimp
            for j = i:dimp
                coefmon = feval(symengine,'poly2list',FPexpr(i,j),charvartable);
                coefmonmatr = subs(coefmon.',decvartable,zeros(1,length(decvartable)));
                for k = 1:size(coefmonmatr,1)
                    s_ijk = coefmonmatr(k,1);

                    dummy_decvar = reshape(coefmon(k),2,1);
                    s_ijk_decvar = dummy_decvar;

                    mon_k= double(coefmonmatr(k,2:end));
                    [val,ind_k] = max(sum((Zfull == kron(ones(nmon,1),mon_k)),2));
                    coeffnts((ind_k-1)*dimp+i,j) = s_ijk;
                    coeffnts((ind_k-1)*dimp+j,i) = s_ijk;

                    coeffnts_decvar((ind_k-1)*dimp+i,j) = s_ijk_decvar;
                    coeffnts_decvar((ind_k-1)*dimp+j,i) = s_ijk_decvar;

                end

            end
        end

    else%JA&GV 07/2013 the code below  addresses the matrix case

        Zfull = [];
        for i = 1:dimp
            for j = i:dimp
                coefmon = feval(symengine,'poly2list',FPexpr(i,j),charvartable);
                coefmon = coefmon.';
                for k = 1:length(coefmon)
                    dummyvar = reshape(coefmon(k),2,1);
                    Z(k,:) = double(dummyvar(2));
                end
                Zfull = sortNoRepeat(Zfull,Z);
                Z = [];
            end
        end
        Z = sparse(Zfull);
        [nmon,nvar] = size(Z);
        coeffnts = sparse(nmon*dimp,dimp);

        for i = 1:dimp
            for j = i:dimp
                coefmon = feval(symengine,'poly2list',FPexpr(i,j),charvartable);
                coefmon = coefmon.';
                for k = 1:length(coefmon)
                    dummyvar = reshape(coefmon(k),2,1);
                    s_ijk = double(dummyvar(1));
                    mon_k= double(dummyvar(2));
                    [val,ind_k] = max(sum((Zfull == kron(ones(nmon,1),mon_k))'));
                    coeffnts((ind_k-1)*dimp+i,j) = s_ijk;
                    coeffnts((ind_k-1)*dimp+j,i) = s_ijk;
                end
            end
        end

    end

    % Constructing At and b matrices
    At = sparse(decvarnum,size(Z,1)*dimp^2);

    %JA&GV may/2013 uses the jacobian function to derive the expression
    if isempty(decvartable)==0
        for i = 1:size(Z,1)
            Mivec = reshape(coeffnts_decvar((i-1)*dimp+1:i*dimp,:),dimp^2,1);
            At(:,(i-1)*dimp^2+1:i*dimp^2) = sparse(-double(jacobian(Mivec,decvartable))');
        end
    end

    b = sparse(coeffnts);
end

author: yu-jiang

repo: radpbook

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


text:
function Score = TwoFoldCV(hfig,numK2,M_stim,timelists_names) % need to LoadFullFish first
%% partitions for CV
% timelists = getappdata(hfig,'timelists');
timelists_names = getappdata(hfig,'timelists_names');
% periods = getappdata(hfig,'periods');
% if length(periods)>1,
%     timelistsCV = cell(length(M_stim),2);
%     for k_stim = 1:length(M_stim), % :3
%         i_stim = M_stim(k_stim);
%         TL = timelists{i_stim};
%         period = periods(i_stim);
%         nrep = size(TL,2)/periods(i_stim); % integer
%         n = floor(nrep/2);
%         timelistsCV{k_stim,1} = TL(1):TL(n*period);
%         timelistsCV{k_stim,2} = TL(1+n*period):TL(2*n*period);
%     end
% end

%% CV loop: auto-clustering with the partitions
absIX = getappdata(hfig,'absIX');

Score = zeros(length(M_stim),2);

for k_stim = 1:length(M_stim), % :3
    i_stim = M_stim(k_stim);
    CIX = cell(1,2);
    GIX = cell(1,2);
    for k = 1:2,
        %% Cluster to start auto-clustering
        i_ClusGroup = M_ClusGroup(k_fish);
        i_Cluster = M_Cluster(k_fish);
        ClusGroup = VAR(i_fish).ClusGroup{i_ClusGroup};
        numK = ClusGroup(i_Cluster).numK;
        gIX = ClusGroup(i_Cluster).gIX;
        cIX_abs = ClusGroup(i_Cluster).cIX_abs; % convert absolute index to index used for this dataset
        [~,cIX] = ismember(cIX_abs,absIX);

        % ~UpdateTimeIndex
        tIX = timelistsCV{k_stim,k};
        M_0 = GetTimeIndexedData_Default_Direct(hfig,cIX,tIX,'isAllCells');

        isWkmeans = 1;
        [cIX,gIX] = AutoClustering(cIX,gIX,absIX,i_fish,M_0,isWkmeans,numK2);
        CIX{k} = cIX;
        GIX{k} = gIX;
    end
    % plot cell-matching figure
    Score(k_stim,1) = HungarianCV(CIX{1},CIX{2},GIX{1},GIX{2},timelists_names{i_stim});
    Score(k_stim,2) = HungarianCV(CIX{2},CIX{1},GIX{2},GIX{1},timelists_names{i_stim});
end
end

author: xiuyechen

repo: FishExplorer


text:
function P = adjustPsfSize(P)
%
%      P = adjustPsfSize(P);
%
%  If the size of a PSF is small compared to the dimensions of
%  the image, then our matrix vector multiply functions will be
%  inefficient due to our overlap add and overlap save approach.
%  Therefore, we make the minimum PSF size to be 32.
%
%  Input:  P - psf object
%  Output: P - psf object
%

%  J. Nagy 1/12/02

I = P.image;
center = P.center;

for k = 1:size(I,3)
  for j = 1:size(I,2)
    for i = 1:size(I,1)
      if length(I{i,j,k}) < 32
        I1 = zeros(32*ones(1,length(size(I{i,j,k}))));
        I1(1:size(I{i,j,k},1), 1:size(I{i,j,k},2), 1:size(I{i,j,k},3)) = I{i,j,k};
        I{i,j,k} = I1;
      end
    end
  end
end
P = psf(I, center);

repo: IRtools


text:
function pass = test_size( pref )
% Test SIZE

if ( nargin == 0 )
    pref = chebfunpref;
end

% Check two components:
F = chebfun2v(@(x,y) cos(x), @(x,y) sin(y));
pass(1) = all( size(F) == [2 inf inf]);
pass(2) = ( size(F, 1) == 2 );
pass(3) = ( size(F, 2) == inf);
pass(4) = ( size(F, 3) == inf);

% Check three components:
F = chebfun2v(@(x,y) cos(x), @(x,y) sin(y), @(x,y) cos(x));
pass(5) = all( size(F) == [3 inf inf]);
pass(6) = ( size(F, 1) == 3 );
pass(7) = ( size(F, 2) == inf);
pass(8) = ( size(F, 3) == inf);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% StackExchange Mathematics Q3537129
% https://math.stackexchange.com/questions/3537129
% Projection of z onto the Affine Half Space {x∣Ax=b,x≥0}.
% References:
%   1.
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     27/12/2020
%   *   First release.

%% General Parameters

subStreamNumberDefault = 2088;0; %<! Set to 0 for Random

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Parameters

numRows         = 10;
numCols         = 20;

% Solvers parameters
numIterations   = 500;
stopThr         = 1e-6;

%% Load / Generate Data

mA = randn(numRows, numCols);
vB = randn(numRows, 1);
vZ = randn(numCols, 1);

mAA     = mA * mA.';
mInvAA  = pinv(mAA);

hObjFun = @(vX) 0.5 * sum((vX - vZ) .^ 2);

hP1 = @(vX) vX - (mA.' * mInvAA * ((mA * vX) - vB)); %<! Projection onto the 1st constraint set
hP2 = @(vX) max(vX, 0); %<! Projection onto the 2nd constraint set

hP = @(vX) OrthogonalProjectionOntoConvexSets({hP1; hP2}, vX, numIterations, stopThr);

solverIdx       = 0;
cMethodString   = {};

%% Solution by CVX

solverString = 'CVX';

% cvx_solver('SDPT3'); %<! Default, Keep numRows low
% cvx_solver('SeDuMi');
% cvx_solver('Mosek'); %<! Can handle numRows > 500, Very Good!
% cvx_solver('Gurobi');

hRunTime = tic();

cvx_begin('quiet')
% cvx_begin()
    % cvx_precision('best');
    variable vX(numCols, 1);
    minimize( 0.5 * sum_square(vX - vZ) );
    subject to
        mA * vX == vB;
        vX >= 0;
cvx_end

runTime = toc(hRunTime);

% vX = mX(:);

disp([' ']);
disp([solverString, ' Solution Summary']);
disp(['The ', solverString, ' Solver Status - ', cvx_status]);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(vX))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX(:).'), ' ]']);
disp(['The Run Time Is Given By - ', num2str(runTime), ' [Sec]']);
disp([' ']);

sCvxSol.vXCvx     = vX;
sCvxSol.cvxOptVal = hObjFun(vX);

%% Solution by Dykstra's Projection Algorithm

solverIdx                   = solverIdx + 1;
cLegendString{solverIdx}    = ['Projected Gradient Method'];

hRunTime = tic();

vX = hP(vZ);

runTime = toc(hRunTime);

disp([' ']);
disp([cLegendString{solverIdx}, ' Solution Summary']);
disp(['The Optimal Value Is Given By - ', num2str(hObjFun(vX))]);
disp(['The Optimal Argument Is Given By - [ ', num2str(vX(:).'), ' ]']);
disp(['The Run Time Is Given By - ', num2str(runTime), ' [Sec]']);
disp([' ']);

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function [newModel, modelChanges] = addMetabolicNetwork(model,rxnsToAdd,metsToAdd)
%addMetabolicNetwork  Integrate a new metabolic network into a model
%
% NOTE: The function is to incorporate a subnetwork that is provided by two
%       arrays, rxnsToAdd and metsToAdd, which should be in defined
%       structure as below，and with reactions and metabolites that are
%       NOT included in the model
%
% Usage:
%   [newModel, modelChanges] = addMetabolicNetwork(model,rxnsToAdd,metsToAdd)
%
%
% Inputs:
%   model        A model structure
%
%   rxnsToAdd    the reaction structure can have the following fields:
%            rxns               cell array with unique strings that
%                               identifies each reaction
%            equations          cell array with equation strings. Decimal
%                               coefficients are expressed as "1.2".
%                               Reversibility is indicated by "<=>" or "=>"
%                               The metabolites have to be written as
%                               "metNames[comps]". Only compartments in
%                               model.comps are allowed. New metabolites to
%                               the model should be added to metsToAdd array
%            rxnNames           cell array with the names of each reaction
%                               (opt, default '')
%            lb                 vector with the lower bounds (opt, default
%                               to model.annotations.defaultLB or -inf for
%                               reversible reactions and 0 for irreversible
%                               when "equations" is used. When "mets" and
%                               "stoichCoeffs" are used it defaults for all
%                               reactions to model.annotations.defaultLB or
%                               -inf)
%            ub                 vector with the upper bounds (opt, default
%                               to model.annotations.defaultUB or inf)
%            eccodes            cell array with the EC-numbers for each
%                               reactions. Delimit several EC-numbers with
%                               ";" (opt, default '')
%            subSystems         cell array with the subsystems for each
%                               reaction (opt, default '')
%            grRules            cell array with the gene-reaction
%                               relationship for each reaction. For example
%                               "(A and B) or (C)" means that the reaction
%                               could be catalyzed by a complex between
%                               A & B or by C on its own. All the genes
%                               have to be present in model.genes. Add
%                               genes with addGenesRaven before calling
%                               this function if needed (opt, default '')
%            rxnReferences      cell array with reaction references (opt,
%                               default '')
%            rxnConfidenceScores   vector with reaction confidence scores
%                               (opt, default NaN)
%
%   metsToAdd   the metabolite structure can have the following fields:
%            mets              cell array with unique strings that
%                              identifies each metabolite (opt, default IDs
%                              of new metabolites are numbered with the
%                              prefix defined below)
%            metNames          cell array with the names of each
%                              metabolite
%            compartments      cell array with the compartment of each
%                              metabolite. Should match model.comps.
%                              If this is a string rather than a cell
%                              array it is assumed that all mets are in
%                              that compartment
%            inchis            cell array with InChI strings for each
%                              metabolite (opt, default '')
%            metFormulas       cell array with the formulas for each of
%                              the metabolites (opt, default '')
%            metCharges        metabolite charge (opt, default NaN)
%
%
% Outputs:
%   newModel      a new model with the added reactons and metabolites
%
%   modelChanges  a cell structure containing the original and new reaction
%                 and/or metabolite properties. Reaction-related changes
%                 are reported in the modelChanges.rxns field, whereas
%                 metabolite-related changes are reported in the
%                 modelChanges.mets field
%

if nargin<2
    error('Missing input!');
end

if nargin<3
    metsToAdd = [];
end

%% Add new metabolites

if isempty(metsToAdd)
    newModel = model;
else
    % the following metabolite fields are required to be included
    requiredMetFields = {'mets','metNames','metFormulas','metCharges','compartments'};
    if ~all(ismember(requiredMetFields, fieldnames(metsToAdd)))
        error('There is missing fields in the metabolite structure to be added.');
    end

    % verify that none of the metabolites exist in the current model
    if any(ismember(metsToAdd.mets, model.mets))
        error('One or more metabolite IDs to be added already exist in the model.');
    end

    % add metabolites to the model
    newModel = addMets(model,metsToAdd);
end

%% Check new reactions

% the following reaction fields are required to be included
requiredRxnFields = {'rxns','equations','subSystems','grRules'};
if ~all(ismember(requiredRxnFields, fieldnames(rxnsToAdd)))
    error('There is missing fields in the metabolite structure to be added.');
end

% verify that none of the reactions exist in the current model
if any(ismember(rxnsToAdd.rxns, model.rxns))
    error('One or more reactions to be added already exist in the model.');
end

% make sure all equations have metabolites in a format of "metNames[comps]"
if ~all(contains(rxnsToAdd.equations, '[') & contains(rxnsToAdd.equations, ']'))
    error('The metabolites in equations should be in a format of "metNames[comps]".');
end

%% Add new genes

% extract genes from grRules
newGenes = getGenesFromGrRules(rxnsToAdd.grRules);

% remove genes that are already existed in model
if any(ismember(newGenes, newModel.genes))
    fprintf('One or more genes to be added already exist in the model，therefore they will be merged.\n');
    newGenes = setdiff(newGenes, newModel.genes);
end

% append new genes to list of model genes
newModel.genes = [newModel.genes; newGenes];

% add new columns to rxnGeneMat will be updated after the new reactions are added below.
newModel.rxnGeneMat(:, end+1:end+numel(newGenes)) = 0;

%% Integrate new subnetwork to the model
newModel = addRxns(newModel, rxnsToAdd, 3);

% update genes and rxnGeneMat with new grRules
[newModel.genes, newModel.rxnGeneMat] = getGenesFromGrRules(newModel.grRules);

%% Document model changes
modelChanges = docModelChanges(model, newModel);

author: SysBioChalmers

repo: Human-GEM


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

% Copyright 2002-2003 The MathWorks, Inc.

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

% Last Modified by GUIDE v2.5 02-Dec-2005 23:18:34

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

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

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

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

% Update handles structure
guidata(hObject, handles);

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

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

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
x0=[0.85 3.2];
[t,x]=ode45(@LV,[0 10],x0);
axes(handles.axes1);
plot(t,x(:,1),'r');hold on
plot(t,x(:,2),'g');
axis([0,10,0,5]);
Xlabel('temps');
Ylabel('X et Y');
hold off
axes(handles.axes2);
plot(x(:,1),x(:,2),'b');hold on
axis([0.5,1,2.5,3.5]);
Xlabel('X');
Ylabel('Y');
hold off


text:
function Z = loadAudio(paths)
% Syntax:   Z = loadAudio(paths);

% Parse inputs
p = numel(paths);

% Load audio
audio = cell(1,p);
for i = 1:p
    path = paths{i};
    [~,~,ext] = fileparts(path);
    if strcmpi(ext,'.wav')
        % Load from .wav file
        audio{i} = audioread(path,'double');
    elseif strcmpi(ext,'.mat')
        % Load from .mat file
        data = load(path);
        fields = fieldnames(data);
        audio{i} = data.(fields{1});
    else
        % Unsupported extension
        error('Unsupported extension ''%s''',ext);
    end
end

% Determine sample length
nSamples = cellfun(@numel,audio);
n = min(nSamples);

% Construct audio matrix
Z = zeros(p,n);
for i = 1:p
    ni     = nSamples(i);
    gapL   = floor(0.5 * (ni - n));
    gapR   = ceil( 0.5 * (ni - n));
    Z(i,:) = audio{i}((gapL + 1):(ni - gapR));
end

% Normalize audio
Z = normalizeAudio(Z);


text:
% Fig. 5.20   Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
% Script to generate Figure 5.20
n=1;
d=[1 8 32 0];
rlocus(n,d)
axis([-10 6 -6 6])
hold on
x=[ 0 -2 0 -2];
y=[0  2*sqrt(3)  0 -2*sqrt(3) ];
plot(x,y)
r=roots([1 4 16]);
plot(r,'*')
title('Fig. 5.20 Locus for L=1/s(s^2+8s+32)')
z=0:.1:.9;
 wn=1:1:10;
 sgrid(z, wn)
 hold off


text:
function [pg, data, imsegs] = msTestImage(im, imsegs, classifiers, nsegments, normalize, ...
    smaps, spdata, adjlist, edata)
%  [pg, data, imsegs] = msTestImage(im, imsegs, classifiers, nsegments, smaps,
%  spdata, adjlist, edata)
%
% Computes the marginals of the labels for the given input image
% spdata, adjlist, edata are optional inputs
% Note: only first three arguments are required

% if imsegs is a cell, it contains the system command to make the
% segmentation image and the filename of the segmentation image

if iscell(imsegs)
    syscall = imsegs{1};
    outfn = imsegs{2};
    system(syscall);
    imsegs = processSuperpixelImage(outfn);
end

if ~exist('normalize', 'var') || isempty(normalize)
    normalize = 1;
end

labelclassifier = classifiers.labelclassifier;
segclassifier = classifiers.segclassifier;

imdata = mcmcComputeImageData(im, imsegs);

if ~exist('spdata', 'var') || isempty(spdata)
    spdata = mcmcGetSuperpixelData(im, imsegs);
end

if ~exist('smaps', 'var') || isempty(smaps)

    if ~exist('adjlist', 'var') || ~exist('edata', 'var') || isempty(adjlist) || isempty(edata)
        [edata, adjlist] = mcmcGetEdgeData(imsegs, spdata);
    end

    eclassifier = classifiers.eclassifier;
    ecal = classifiers.ecal;
    if isfield(classifiers, 'vclassifierSP')
        vclassifierSP = classifiers.vclassifierSP;
        hclassifierSP = classifiers.hclassifierSP;
        [pvSP, phSP, pE] = mcmcInitialize(spdata, edata, ...
            adjlist, imsegs, vclassifierSP, hclassifierSP, eclassifier, ecal, 'none');
    else
        pE = test_boosted_dt_mc(eclassifier, edata);
        pE = 1 ./ (1+exp(ecal(1)*pE+ecal(2)));
    end
    smaps = generateMultipleSegmentations2(pE, adjlist, imsegs.nseg, nsegments);
end

labdata = cell(1, size(smaps, 2));
for k = 1:size(smaps, 2)
    labdata{1, k} = mcmcGetSegmentFeatures(imsegs, spdata, imdata, smaps(:, k), (1:max(smaps(:, k))));
end

pg = msTest(imsegs, labdata, {smaps}, labelclassifier, segclassifier, normalize);

% nsp = imsegs.nseg;
%
% nclasses = size(labelclassifier.wcs, 2);
% pg = zeros(nsp, nclasses);
%
% segs = cell(nsp, 1);
%
% for k = 1:size(smaps, 2)
%
%     for s = 1:max(smaps(:, k))
%
%         [segs, ind] = checksegs(segs, smaps(:, k), s);
%
%         if ~isempty(ind)
%
%             labdata = mcmcGetSegmentFeatures(imsegs, spdata, imdata, smaps(:, k), s);
%
%             vconf = test_boosted_dt_mc(vclassifier, labdata);
%             vconf = 1 ./ (1+exp(-vconf));
%             vconf = vconf / sum(vconf);
%
%             hconf = test_boosted_dt_mc(hclassifier, labdata);
%             hconf = 1 ./ (1+exp(-hconf));
%             hconf = hconf / sum(hconf);
%
%             sconf = test_boosted_dt_mc(sclassifier, labdata);
%             sconf = 1 ./ (1+exp(-sconf));
%
%             pgs = [vconf(1) vconf(2)*hconf vconf(3)]*sconf;
%
%             pg(ind, :) = pg(ind, :) + repmat(pgs, numel(ind), 1);
%         end
%
%     end
%
% end
%
% pg = pg ./ max(repmat(sum(pg, 2), 1, size(pg, 2)), 0.00001);
pg = pg{1};

data.smaps = smaps;
data.edata = edata;
data.adjlist = adjlist;
data.spdata = spdata;
data.imdata = imdata;

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

function [segs, ind] = checksegs(segs, map, s)
% Checks whether this segment has been seen before

ind = find(map==s);

if isempty(ind)
    return;
end

% if numel(ind)==1 % already accounted for by superpixels
%     ind = [];
%     return;
% end

oldsegs = segs{ind(1)};

for k = 1:numel(oldsegs)
    if (numel(oldsegs{k})==numel(ind)) && all(oldsegs{k}==ind)
        ind = [];
        return;
    end
end

segs{ind(1)}{end+1} = ind;

author: Cloud-CV

repo: object-proposals


text:
% Copyright (C) 2018  Symeon Symeonidis, Stefanos Tsantilas, Stelios Mitilineos
% simos421@gmail.com, steftsantilas@gmail.com, smitil@gmail.com
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.

function CstDefineEfieldMonitor(mws,Efieldname, frequency)

%'@ define monitor: e-field (f=2.25)
% Efieldname is a string

Monitor = invoke(mws,'Monitor');
invoke(Monitor,'Reset');
invoke(Monitor,'Name',Efieldname);
invoke(Monitor,'Dimension','Volume');
invoke(Monitor,'Domain','Frequency');
invoke(Monitor,'FieldType','Efield');
invoke(Monitor,'Frequency',num2str(frequency));
invoke(Monitor,'UseSubvolume','False');
invoke(Monitor,'SetSubvolume','-53.310273111111', '53.310273111111', '-53.310273111111', '53.310273111111', '-33.310273111111', '71.310273111111');
invoke(Monitor,'Create');

end

repo: CST-MATLAB-API


text:
function [ fn ] = eph_name( denum, year )
%eph_name creates a naming convection for the file names of the binary ephemerides

  fn = Ephem.asc_name( denum, year );
  if ~isempty(fn)
    fn(1:3) = 'eph';
  end
end

author: goGPS-Project


text:
% SPM5 UPDATE 23/11/07
% Sets the default values for the FieldMap toolbox
%
% FORMAT pm_defaults_Trio_eFoV_dn
%_______________________________________________________________________
%
% This file is intended for use with the Siemens fieldmap sequence
% on the Trio scanner at the AMRIG/FIL and the eFoV EPI sequences with
% PE blips=+1:
%
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Chloe Hutton and Jesper Andersson
% $Id: pm_defaults_Trio_eFoV_dn.m 5015 2012-10-24 13:40:07Z guillaume $

global pm_def

% Defaults for creating field map. (See pm_make_fieldmap.m and
%                                   FieldMap.man for more info.)
%=======================================================================
pm_def.INPUT_DATA_FORMAT = 'PM';      % 'RI' = load two real and
                                      % imaginary image pairs
                                      % 'PM' = load one or two
                                      % phase and magnitude image
                                      % pairs.
pm_def.SHORT_ECHO_TIME = 10.0;        % Short echo time in ms for Trio
pm_def.LONG_ECHO_TIME = 12.46;        % Long echo time in ms for Trio
pm_def.MASKBRAIN = 1;                 % Do brain masking (1 or 0,
                      % 0 for EPI fieldmaps)

% Defaults for unwrapping options. (See pm_make_fieldmap.m and
%                                   FieldMap.man for more info.)
%=======================================================================
pm_def.UNWRAPPING_METHOD = 'Mark3D';  % Unwrapping options are:
                                      % 'Huttonish', 'Mark3D' or 'Mark2D'
pm_def.FWHM = 10;                     % FWHM of Gaussian filter used to
                                      % implement weighted smoothing of
                                      % unwrapped maps.
pm_def.PAD = 0;                       % Size of padding kernel if required.
pm_def.WS = 1;                        % Weighted or normal smoothing.

% Flags for brain extraction
%=======================================================================
pm_def.MFLAGS.TEMPLATE = fullfile(spm('Dir'),'toolbox','FieldMap','T1.nii');
pm_def.MFLAGS.FWHM = 5;     % In mm
pm_def.MFLAGS.NERODE = 2;   % In voxels
pm_def.MFLAGS.NDILATE = 4;  % In voxels
pm_def.MFLAGS.THRESH = 0.5;
pm_def.MFLAGS.REG = 0.02;   % A larger value helps segmentation to converge
pm_def.MFLAGS.GRAPHICS = 0; % A larger value helps segmentation to converge

% Defaults for converting field map to voxel displacement map.
%=======================================================================
pm_def.EPI_BASED_FIELDMAPS = 0;         % EPI=1, other=0.
pm_def.K_SPACE_TRAVERSAL_BLIP_DIR = +1; % +ve k-space = 1, -ve = -1.
pm_def.TOTAL_EPI_READOUT_TIME = 37.0;   % Trio eFoV EPI RO time

% Defaults for Unwarping.
%=======================================================================
pm_def.DO_JACOBIAN_MODULATION = 0;    % Do jacobian modulation to adjust
                                      % for compression or stretching
                                      % No = 0, Yes = 1


text:
function co = convexification(varargin)
%CONVEXIFICATION Compute the convexification of a polygon.
%
%   CO = convexification(H)
%   Creates convexification from support function. Support function is
%   supposed to be uniformly distributed over [0 2pi].
%
%   CO = convexification(POLYGON)
%   Computes support function of the polygon, then the corresponding
%   convexification.
%
%   CO = convexification(POLYGON, N)
%   Uses N points for convexification computation. Note that the number of
%   points of CO can be lower than N.
%
%   CAUTION: The result will be valid only for convex polygons.
%
%   See also
%   polygons2d, supportFunction

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2005-01-12
% Copyright 2005-2022 INRA - Cepia Software Platform

if ~isempty(varargin)>0
    var = varargin{1};
    if size(var, 2)==1
        h = var;
    else
        poly = var;
        N = 128;
        if length(varargin)>1
            N = varargin{2};
        end
        h = supportFunction(poly, N);
    end
else
    error('not enough input arguments');
end

N   = length(h);
u   = (0:2*pi/N:2*pi*(1-1/N))';
v   = [cos(u) sin(u)].*[h h];

i1  = 1:N;
i2  = [2:N 1];
i3  = [3:N 1 2];

circ = zeros(N, 4);
for i=1:N
    circ(i, 1:4) = createDirectedCircle(v(i1(i),:), v(i2(i),:), v(i3(i), :));
end

% remove non direct-oriented circles
circ = circ(circ(:,4)==0, :);

% keep only circles seen several times
dp = diff(circ(:,1:2));
dp = sum(dp.*dp, 2);
ind1 = [1; find(dp<1e-10)+1];
circ = circ(ind1, :);

% keep only one instance of each circle
dp = diff(circ(:,1:2));
dp = sum(dp.*dp, 2);
ind = [1; find(dp>1e-10)+1];
co = 2*circ(ind, 1:2);

% eventually remove the last point if it is the same as the first one
if distancePoints(co(1,:), co(end, :))<1e-10 && size(co, 1)>1
    co = co(1:end-1,:);
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% loadtxt() - load ascii text file into numeric or cell arrays
%
% Usage:
%   >> array = loadtxt( filename, 'key', 'val' ...);
%
% Inputs:
%    filename - name of the input file
%
% Optional inputs
%   'skipline' - number of lines to skip {default:0}. If this number is
%                negative the program will only skip non-empty lines
%                (can be usefull for files transmitted from one platform
%                to an other, as CR may be inserted at every lines).
%   'convert'  - 'on' standard text conversion, see note 1
%                'off' no conversion, considers text only
%                'force' force conversion, NaN are returned
%                for non-numeric inputs {default:'on'}
%   'delim'    - ascii character for delimiters. {default:[9 32]
%                i.e space and tab}. It is also possible to enter
%                strings, Ex: [9 ' ' ','].
%   'blankcell' - ['on'|'off'] extract blank cells {default:'on'}
%   'verbose'  - ['on'|'off'] {default:'on'}
%   'convertmethod' - ['str2double'|'str2num'] default is 'str2double'
%   'nlines'   - [integer] number of lines to read {default: all file}
%
% Outputs:
%    array - cell array. If the option 'force' is given, the function
%            retrun a numeric array.
%
% Notes: 1) Since it uses cell arrays, the function can handle text input.
%        The function reads each token and then try to convert it to a
%        number. If the conversion is unsucessfull, the string itself
%        is included in the array.
%        2) The function adds empty entries for rows that contains
%        fewer columns than others.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 29 March 2002

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

function array = loadtxt( filename, varargin );

if nargin < 1
	help loadtxt;
	return;
end;
if ~isempty(varargin)
   try, g = struct(varargin{:});
   catch, disp('Wrong syntax in function arguments'); return; end;
else
    g = [];
end;

g = finputcheck( varargin, { 'convert'   'string'   { 'on';'off';'force' }   'on';
                             'skipline'  'integer'  [0 Inf]          0;
                             'verbose'   'string'   { 'on';'off' }   'on';
                             'uniformdelim' 'string'   { 'on';'off' }   'off';
                             'blankcell' 'string'   { 'on';'off' }   'on';
                             'convertmethod' 'string'   { 'str2double';'str2num' }   'str2double';
                             'delim'     { 'integer';'string' } []               [9 32];
                             'nlines'    'integer'  []               Inf });
if isstr(g), error(g); end;
if strcmpi(g.blankcell, 'off'), g.uniformdelim = 'on'; end;
g.convert = lower(g.convert);
g.verbose = lower(g.verbose);
g.delim = char(g.delim);

% open the file
% -------------
if exist(filename) ~=2, error( ['file ' filename ' not found'] ); end;
fid=fopen(filename,'r','ieee-le');
if fid<0, error( ['file ' filename ' found but error while opening file'] ); end;

index = 0;
while index < abs(g.skipline)
    tmpline = fgetl(fid);
    if g.skipline > 0 | ~isempty(tmpline)
        index = index + 1;
    end;
end; % skip lines ---------

inputline = fgetl(fid);
linenb = 1;
if strcmp(g.verbose, 'on'), fprintf('Reading file (lines): '); end;
while isempty(inputline) | inputline~=-1
     colnb = 1;
     if ~isempty(inputline)
         tabFirstpos = 1;

         % convert all delimiter to the first one
         if strcmpi(g.uniformdelim, 'on')
             for index = 2:length(g.delim)
                 inputline(find(inputline == g.delim(index))) = g.delim(1);
             end;
         end;

         while ~isempty(deblank(inputline))
             if strcmpi(g.blankcell,'off'), inputline = strtrim(inputline); end;
             if tabFirstpos && length(inputline) > 1 && all(inputline(1) ~= g.delim), tabFirstpos = 0; end;
             [tmp inputline tabFirstpos] = mystrtok(inputline, g.delim, tabFirstpos);
             switch g.convert
                case 'off', array{linenb, colnb} = tmp;
                case 'on',
                     if strcmpi(g.convertmethod, 'str2double')
                         tmp2 = str2double(tmp);
                         if isnan( tmp2 )  , array{linenb, colnb} = tmp;
                         else                array{linenb, colnb} = tmp2;
                         end;
                     else
                         tmp2 = str2num(tmp);
                         if isempty( tmp2 )  , array{linenb, colnb} = tmp;
                         else                  array{linenb, colnb} = tmp2;
                         end;
                     end;
                case 'force', array{linenb, colnb} = str2double(tmp);
             end;
             colnb = colnb+1;
         end;
	     linenb = linenb +1;
     end;
     inputline = fgetl(fid);
     if linenb > g.nlines
         inputline = -1;
     end;
     if ~mod(linenb,10) & strcmp(g.verbose, 'on'), fprintf('%d ', linenb); end;
end;
if strcmp(g.verbose, 'on'),  fprintf('%d\n', linenb-1); end;
if strcmp(g.convert, 'force'), array = [ array{:} ]; end;
fclose(fid);

% problem strtok do not consider tabulation
% -----------------------------------------
function [str, strout, tabFirstpos] = mystrtok(strin, delim, tabFirstpos);
    % remove extra spaces at the beginning
    while any(strin(1) == delim) && strin(1) ~= 9 && strin(1) ~= ','
         strin = strin(2:end);
    end;
    % for tab and coma, consider empty cells
    if length(strin) > 1 && any(strin(1) == delim)
        if tabFirstpos || any(strin(2) == delim)
            str = '';
            strout = strin(2:end);
            if strin(2) ~= 9 && strin(2) ~= ','
                tabFirstpos = 0;
                strout = strtrim(strout);
            end;
        else
            [str, strout] = strtok(strin, delim);
        end;
    else
        [str, strout] = strtok(strin, delim);
    end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function prob_test054 ( )

%*****************************************************************************80
%
%% TEST054 tests EMPIRICAL_DISCRETE_CDF, EMPIRICAL_DISCRETE_CDF_INV, EMPIRICAL_DISCRETE_PDF;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  a = 6;

  b(1:a) = [ 1.0, 1.0, 3.0, 2.0, 1.0, 2.0 ];
  c(1:a) = [ 0.0, 1.0, 2.0, 4.5, 6.0, 10.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST054\n' );
  fprintf ( 1, '  For the Empirical Discrete PDF:\n' );
  fprintf ( 1, '  EMPIRICAL_DISCRETE_CDF evaluates the CDF;\n' );
  fprintf ( 1, '  EMPIRICAL_DISCRETE_CDF_INV inverts the CDF.\n' );
  fprintf ( 1, '  EMPIRICAL_DISCRETE_PDF evaluates the PDF;\n' );

  check = empirical_discrete_check ( a, b, c );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST054 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A = %6d\n', a );
  r8vec_print ( a, b, '  PDF parameter B:' );
  r8vec_print ( a, c, '  PDF parameter C:' );

  seed = 123456789;

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

  for i = 1 : 10

    [ x, seed ] = empirical_discrete_sample ( a, b, c, seed );

    pdf = empirical_discrete_pdf ( x, a, b, c );

    cdf = empirical_discrete_cdf ( x, a, b, c );

    x2 = empirical_discrete_cdf_inv ( cdf, a, b, c );

    fprintf ( 1, ' %14f  %14f  %14f  %14f\n', x, pdf, cdf, x2 );

  end

  return
end


text:
function wMed = weightedMedian(D,W)

% ----------------------------------------------------------------------
% Function for calculating the weighted median
% Sven Haase
%
% For n numbers x_1,...,x_n with positive weights w_1,...,w_n,
% (sum of all weights equal to one) the weighted median is defined as
% the element x_k, such that:
%           --                        --
%           )   w_i  <= 1/2   and     )   w_i <= 1/2
%           --                        --
%        x_i < x_k                 x_i > x_k
%
%
% Input:    D ... matrix of observed values
%           W ... matrix of weights, W = ( w_ij )
% Output:   wMed ... weighted median
% ----------------------------------------------------------------------

if nargin ~= 2
    error('weightedMedian:wrongNumberOfArguments', ...
      'Wrong number of arguments.');
end

if size(D) ~= size(W)
    error('weightedMedian:wrongMatrixDimension', ...
      'The dimensions of the input-matrices must match.');
end

% normalize the weights, such that: sum ( w_ij ) = 1
% (sum of all weights equal to one)

WSum = sum(W(:));
W = W / WSum;

% (line by line) transformation of the input-matrices to line-vectors
d = reshape(D',1,[]);
w = reshape(W',1,[]);

% sort the vectors
A = [d' w'];
ASort = sortrows(A,1);

dSort = ASort(:,1)';
wSort = ASort(:,2)';
sumVec = cumsum(wSort);

wMed = [];
j = 0;

while isempty(wMed)
    j = j + 1;
    if sumVec(j) >= 0.5
        wMed = dSort(j);    % value of the weighted median
    end
end

% final test to exclude errors in calculation
if ( sum(wSort(1:j-1)) > 0.5 ) & ( sum(wSort(j+1:length(wSort))) > 0.5 )
     error('weightedMedian:unknownError', ...
      'The weighted median could not be calculated.');
end


text:
function x = uncertain(x,varargin)
%UNCERTAIN Declares a variable as uncertain
%
%   F = UNCERTAIN(W) is used to describe the set of uncertain variables
%   in an uncertain program
%
%   INPUT
%    W : SDPVAR object or list of constraints
%    S : Optional distribution information for random uncertainty
%
%   OUTPUT
%    F : Constraint object
%
%   Uncertain is used to declare uncertain variables in robust
%   deterministic worst-case optimization. It can also be used to specify
%   uncertain random variables, and their associated distribution, to be
%   used in OPTIMIZER objects with the SAMPLE command.
%
%   EXAMPLE
%
%    Robust worst-case optimization
%
%    sdpvar x w
%    F = [x + w <= 1], W = [-0.5 <= w <= 0.5];
%    optimize([F,W,uncertain(w)],-x)
%
%    sdpvar x w
%    F = [x + w <= 1], W = [-0.5 <= w <= 0.5];
%    optimize([F,uncertain(W)],-x)
%
%   To specify random uncertainties, you specify the distribution, and all
%   distribution parameters following the syntax n the RANDOM command in
%   the Statistics Toolbox
%
%    sdpvar x w
%    F = [x + w <= 1, uncertain(w, 'uniform',0,1)];
%    P = optimizer([F,W,uncertain(w)],-x,[],w,x)
%    S = sample(P,10); % Sample ten instances and concatenate models
%    S([])             % Solve and return optimal x
%
%    Alternatively, you can specify a function handle which generates
%    samples. YALMIP will always send a trailing argument with dimensions
%
%    F = [x + w <= 1, uncertain(w,@mysampler,myarguments1,...)];
%
%    The standard random case above would thus be recovered with
%
%    F = [x + w <= 1, uncertain(w,@random,0,1)];
%
%
%   See also OPTIMIZE, ROBUSTMODEL, OPTIMIZER, SAMPLE

if nargin == 1 || ((nargin == 2) && strcmpi(varargin{1},'deterministic'))
    x.typeflag = 15;
    x.extra.distribution.name = 'deterministic';
    x = lmi(x);
else
    x.typeflag = 16;
    if isa(varargin{1},'function_handle')
         temp = {varargin{:},x.dim};
    else
        temp = {@random,varargin{:},x.dim};
    end
    x.extra.distribution.name = temp{1};
    x.extra.distribution.parameters = {temp{2:end-1}};
    try
        temp = feval(temp{:});
    catch
        disp(lasterr);
        error('Trial evaluation of attached sample generator failed.')
    end
    x = lmi(x);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function a = read_cor(fname)
% a = read_cor(fname)

%
% read_cor.m
%
% Original Author: Bruce Fischl
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:12 $
%    $Revision: 1.3 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

a = zeros(256, 256, 256);

for i=1:256
%	disp(i);
	fname1 = sprintf('%s/COR-%03d', fname, i);
	fid = fopen(fname1, 'r');
	a(:, :, i) = fread(fid, [256 256], 'uchar');
	fclose(fid);
end

% eof

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function varargout = sl_synch_gui(varargin)
% SL_SYNCH_GUI M-file for sl_synch_gui.fig
%      SL_SYNCH_GUI, by itself, creates a new SL_SYNCH_GUI or raises the existing
%      singleton*.
%
%      H = SL_SYNCH_GUI returns the handle to a new SL_SYNCH_GUI or the handle to
%      the existing singleton*.
%
%      SL_SYNCH_GUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in SL_SYNCH_GUI.M with the given input arguments.
%
%      SL_SYNCH_GUI('Property','Value',...) creates a new SL_SYNCH_GUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before sl_synch_gui_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to sl_synch_gui_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES
%
% Edited by Will Campbell, MathWorks Inc.
% will.campbell@mathworks.com
% January 7, 2010
% Copyright 2010 The MathWorks, Inc.

% Last Modified by GUIDE v2.5 19-Feb-2010 10:51:33

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

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

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

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

% Update handles structure
guidata(hObject, handles);

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

% Get the value of the manual switch and synch the GUI
value = get_param([bdroot '/Manual Switch'],'sw');
if strcmp(value,'0')
    set(handles.checkbox_switch,'Value',0)
else
    set(handles.checkbox_switch,'Value',1)
end

% Get the value of the gain edit box and synch the GUI
value = get_param([bdroot '/Gain'],'Gain');
set(handles.edit_Gain,'String',value);

% Set the value of the gain slider, with max/min of +10/-10
slider_position = max(0,min(1,(str2double(value) + 10)/20));
set(handles.slider_Gain,'Value',slider_position);

% Assign gui, startstop, switch, and gain handles to the base workspace
assignin('base','sl_synch_handles',handles)
assignin('base','startstop_hObject',handles.pushbutton_startstop)
assignin('base','switch_hObject',handles.checkbox_switch)
assignin('base','gain_hObject',handles.edit_Gain)

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

% Get default command line output from handles structure
varargout{1} = handles.output;

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

mystring = get(hObject,'String');
status = get_param(bdroot,'simulationstatus');

if strcmp(mystring,'Start Simulation')

    % Check the status of the simulation and start it if it's stopped
    if strcmp(status,'stopped')
        set_param(bdroot,'simulationcommand','start')
    end

    % Update the string on the pushbutton
    set(handles.pushbutton_startstop,'String','Stop Simulation')

elseif strcmp(mystring,'Stop Simulation')

    % Check the status of the simulation and stop it if it's running
    if strcmp(status,'running')
        set_param(bdroot, 'SimulationCommand', 'Stop')
    end

    % Update the string on the pushbutton
    set(handles.pushbutton_startstop,'String','Start Simulation')

else
    warning('Unrecognized string for pushbutton_startstop') %#ok<WNTAG>
end

% Assign handles and the startstop object to the base workspace
assignin('base','sl_synch_handles',handles)
assignin('base','startstop_hObject',handles.pushbutton_startstop)

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

% Hint: get(hObject,'Value') returns toggle state of checkbox_switch

value = get(hObject,'Value');

if value == 0
    set_param([bdroot '/Manual Switch'],'sw','0')
else
    set_param([bdroot '/Manual Switch'],'sw','1')
end

% Assign handles and the switch object to the base workspace
assignin('base','sl_synch_handles',handles)
assignin('base','switch_hObject',handles.checkbox_switch)

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

% Hints: get(hObject,'String') returns contents of edit_Gain as text
%        str2double(get(hObject,'String')) returns contents of edit_Gain as a double

value = get(hObject,'String');

% Update the model's gain value
set_param([bdroot '/Gain'],'Gain',value)

% Set the value of the gain slider, with max/min of +10/-10
slider_position = max(0,min(1,(str2double(value) + 10)/20));
set(handles.slider_Gain,'Value',slider_position);

% Update simulation if the model is running
status = get_param(bdroot,'simulationstatus');
if strcmp(status,'running')
    set_param(bdroot, 'SimulationCommand', 'Update')
end

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

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

% --- Executes on slider movement.
function slider_Gain_Callback(hObject, eventdata, handles)
% hObject    handle to slider_Gain (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

% Get the value of the gain slider and determine what the gain value should be
slider_position = get(hObject,'Value');
value = num2str(slider_position*20 - 10);

% Update the model's gain value
set_param([bdroot '/Gain'],'Gain',value)

% Set the value of the gain edit box
set(handles.edit_Gain,'String',value);

% Update simulation if the model is running
status = get_param(bdroot,'simulationstatus');
if strcmp(status,'running')
    set_param(bdroot, 'SimulationCommand', 'Update')
end

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

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

author: Sable

repo: mcbench-benchmarks


text:
function R=return_Rt_matrix(alpha,beta,gamma,tx,ty,tz)

% Copyright (C) <2007>  <Francesc Moreno-Noguer, Vincent Lepetit, Pascal Fua>
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the version 3 of the GNU General Public License
% as published by the Free Software Foundation.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
%
% Francesc Moreno-Noguer, CVLab-EPFL, September 2007.
% fmorenoguer@gmail.com, http://cvlab.epfl.ch/~fmoreno/

R(1,1)=cos(alpha)*cos(gamma)-cos(beta)*sin(alpha)*sin(gamma);
R(2,1)=cos(gamma)*sin(alpha)+cos(alpha)*cos(beta)*sin(gamma);
R(3,1)=sin(beta)*sin(gamma);
R(4,1)=0;

R(1,2)=-cos(beta)*cos(gamma)*sin(alpha)-cos(alpha)*sin(gamma);
R(2,2)=cos(alpha)*cos(beta)*cos(gamma)-sin(alpha)*sin(gamma);
R(3,2)=cos(gamma)*sin(beta);
R(4,2)=0;

R(1,3)=sin(alpha)*sin(beta);
R(2,3)=-cos(alpha)*sin(beta);
R(3,3)=cos(beta);
R(4,3)=0;

R(:,4)=[tx,ty,tz,1]';

author: cvlab-epfl

repo: EPnP

save_path: github-repos/MATLAB/cvlab-epfl-EPnP


text:
function [p,t] = smoothmesh(p,t,maxit,tol)

%*****************************************************************************80
%
%  SMOOTHMESH: Smooth a triangular mesh using Laplacian smoothing.
%
% Laplacian smoothing is an iterative process that generally leads to an
% improvement in the quality of the elements in a triangular mesh.
%
%  [p,t] = smoothmesh(p,t);
%
%  p     : Nx2 array of nodal XY coordinates, [x1,y1; x2,y2; etc].
%  t     : Mx3 array of triangles as indices, [n11,n12,n13;
%                                              n21,n22,n23; etc].
%  maxit : Maximum allowable iterations.
%  tol   : Convergence tolerance (Percentage change in edge length must be
%          less than TOL).
%
% If MAXIT or TOL are left empty the default values MAXIT = 20 and TOL =
% 0.01 are used.
%
%  EXAMPLE:
%
%  [p,t] = smoothmesh(p,t,10,0.05);
%
%  Author:
%
%    Darren Engwirda
%

if nargin<4
   tol = [];
   if nargin<3
      maxit = [];
      if nargin<2
         error('Incorrect number of inputs.');
      end
   end
elseif nargin>5
   error('Incorrect number of inputs.')
end
if nargout>2
   error('Incorrect number of outputs.');
end
if isempty(tol)
   tol = 0.01;
end
if isempty(maxit)
   maxit = 20;
end

[p,t] = fixmesh(p,t);                                                      % Ensure consistent mesh

n = size(p,1);
S = sparse(t(:,[1,1,2,2,3,3]),t(:,[2,3,1,3,1,2]),1,n,n);                   % Sparse connectiity matrix
W = sum(S,2);
if any(W==0)
   error('Invalid mesh. Hanging nodes found.');
end

edge = [t(:,[1,2]); t(:,[1,3]); t(:,[2,3])];                               % Non-unique edges
edge = sortrows(sort(edge,2));                                             % Put shared edges next to each other
idx = all(diff(edge)==0,2);                                                % Find shared edges
idx = [idx;false]|[false;idx];                                             % True for all shared edges
bnde = edge(~idx,:);                                                       % Boundary edges
edge = edge(idx,:);                                                        % Internal edges
edge = [bnde; edge(1:2:end-1,:)];                                          % Unique edges
bnd = unique(bnde);                                                        % Boundary nodes

L = max(sqrt(sum((p(edge(:,1),:)-p(edge(:,2),:)).^2,2)),eps);              % Edge length

for iter = 1:maxit
   pnew = (S*p)./[W,W];                                                    % Laplacian smoothing
   pnew(bnd,:) = p(bnd,:);                                                 % Dont let BND nodes move
   p = pnew;

   Lnew = max(sqrt(sum((p(edge(:,1),:)-p(edge(:,2),:)).^2,2)),eps);        % Edge length
   move = norm((Lnew-L)./Lnew,inf);                                        % Percentage change
   if move<tol
      break
   end
   L = Lnew;
end
if iter==maxit
   disp('WARNING: Maximum number of iterations reached, solution did not converge!');
end

end      % smoothmesh()

author: johannesgerer

repo: jburkardt-m


text:
function [assignment, cost] = assignmentsuboptimal1(distMatrix)
%ASSIGNMENTSUBOPTIMAL1    Compute suboptimal assignment
%   ASSIGNMENTSUBOPTIMAL1(DISTMATRIX) computes a suboptimal assignment
%   (minimum overall costs) for the given rectangular distance or cost
%   matrix, for example the assignment of tracks (in rows) to observations
%   (in columns). The result is a column vector containing the assigned
%   column number in each row (or 0 if no assignment could be done).
%
%   [ASSIGNMENT, COST] = ASSIGNMENTSUBOPTIMAL1(DISTMATRIX) returns the
%   assignment vector and the overall cost.
%
%   The algorithm is designed for distance matrices with many forbidden and
%   singly validated assignments (rows or columns containing only one
%   finite element). The algorithm first searches the matrix for singly
%   validated columns and rejects all assignments with multiply validated
%   rows. Afterwards, singly validated rows are searched and assignments to
%   multiply validated columns are rejected. Then, for each row that
%   validates only with singly validated columns (and the other way
%   around), the minimum element is chosen and the assignment is made. If
%   there are still assignments open, the minimum element in the distance
%   matrix is searched and the corresponding assignment is made.
%
%   In scenarios without any forbidden assignments, the algorithm reduces
%   to the last step, which will provide the same result as ASSIGNMENTOPTIMAL2.
%   If there are only some assignments forbidden, the algorithm will perform
%   poorly because singly validated assignments are preferred.
%
%   The last step can still be optimized, see the comments in
%   ASSIGNMENTOPTIMAL2.
%
%   assignment.html  File Exchange  Donate via PayPal
%
%   Markus Buehren
%   Last modified 05.07.2011

% initialize
[nOfRows, nOfColumns] = size(distMatrix);
nOfValidObservations  = zeros(nOfRows,1);
nOfValidTracks        = zeros(1,nOfColumns);
assignment            = zeros(nOfRows,1);
cost                  = 0;

% compute number of validations for each track
for row=1:nOfRows
  nOfValidObservations(row) = length(find(isfinite(distMatrix(row,:))));
end

if any(nOfValidObservations < nOfColumns)

  if all(nOfValidObservations == 0)
    return
  end

  repeatSteps = 1;
  while repeatSteps

    repeatSteps = 0;

    % step 1: reject assignments of multiply validated tracks to singly validated observations
    for col=1:nOfColumns
      index = isfinite(distMatrix(:,col));
      nOfValidTracks(col) = length(find(index));
      if any(nOfValidObservations(index) == 1)
        index = index & (nOfValidObservations > 1);
        if any(index)
          distMatrix(index, col)      = inf;
          nOfValidObservations(index) = nOfValidObservations(index) - 1;
          nOfValidTracks(col)         = nOfValidTracks(col) - length(find(index));
          repeatSteps = 1;
        end
      end
    end

    % step 2: reject assignments of multiply validated observations to singly validated tracks
    if nOfColumns > 1
      for row=1:nOfRows
        index = isfinite(distMatrix(row,:));
        if any(nOfValidTracks(index) == 1)
          index = index & (nOfValidTracks > 1);
          if any(index)
            distMatrix(row, index)    = inf;
            nOfValidTracks(index)     = nOfValidTracks(index) - 1;
            nOfValidObservations(row) = nOfValidObservations(row) - length(find(index));
            repeatSteps = 1;
          end
        end
      end
    end

  end % while repeatSteps
      %disp(sprintf('xx = %d', xx));

  % for each multiply validated track that validates only with singly validated
  % observations, choose the observation with minimum distance
  for row=1:nOfRows
    if nOfValidObservations(row) > 1
      index = isfinite(distMatrix(row,:));
      if all(nOfValidTracks(index) == 1)
        [minDist, col] = min(distMatrix(row,:));
        assignment(row)    = col;
        cost               = cost + minDist;
        distMatrix(row,:)  = inf;
        distMatrix(:,col)  = inf;
      end
    end
  end

  % for each multiply validated observation that validates only with singly validated
  % tracks, choose the track with minimum distance
  for col=1:nOfColumns
    if nOfValidTracks(col) > 1
      index = isfinite(distMatrix(:,col));
      if all(nOfValidObservations(index) == 1)
        [minDist, row] = min(distMatrix(:,col));
        assignment(row)    = col;
        cost               = cost + minDist;
        distMatrix(row,:)  = inf;
        distMatrix(:,col)  = inf;
      end
    end
  end

end

% now, recursively search for the minimum element and do the assignment
while 1

  % find minimum distance observation-to-track pair
  [minDist, index1] = min(distMatrix, [], 1);
  [minDist, index2] = min(minDist);
  row = index1(index2);
  col = index2;

  if isfinite(minDist)

    % make the assignment
    assignment(row)    = col;
    cost               = cost + minDist;
    distMatrix(row, :) = inf;
    distMatrix(:, col) = inf;

  else
    break
  end

end

author: GenLouvain

repo: GenLouvain


text:
% Fig. 5.2   Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
clf
hold off
n2=1;
d2=[1 1 0];
rlocus(n2,d2)
axis([-2 2 -1.5 1.5])
title('Fig.5.2 Root locus of 1/s(s+1)')
grid on
  z=0:.1:.9;
 wn=.5:.5:2;
 sgrid(z, wn)


text:
function W = get_AWGN_transition_probability(sigma, v)
alpha = get_Clambda_zero_points(v);
y = get_y_interval(sigma, alpha);
W = upgrading_transform_AWGN_to_DMC(y, alpha, sigma, v);
end

author: YuYongRun

repo: PolarCodeDecodersInMatlab