text:
function diftab = dif_basis_i ( ival, ntab, xtab )

%*****************************************************************************80
%
%% DIF_BASIS_I computes the I-th Lagrange basis polynomial in divided difference form.
%
%  Discussion:
%
%    The I-th Lagrange basis polynomial for a set of NTAB X values XTAB,
%    L(I,NTAB,XTAB)(X) is a polynomial of order NTAB-1 which is zero at
%    XTAB(J) for J not equal to I, and 1 when J is equal to I.
%
%    The Lagrange basis polynomials have the property that the interpolating
%    polynomial through a set of NTAB data points (XTAB,YTAB) may be
%    represented as
%
%      P(X) = Sum ( 1 <= I <= N ) YTAB(I) * L(I,NTAB,XTAB)(X)
%
%    Higher order interpolation at selected points may be accomplished
%    using repeated X values, and scaled derivative values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Carl deBoor,
%    A Practical Guide to Splines,
%    Springer, 2001,
%    ISBN: 0387953663,
%    LC: QA1.A647.v27.
%
%  Parameters:
%
%    Input, integer IVAL, the index of the desired Lagrange basis polynomial.
%    IVAL should be between 1 and NTAB.
%
%    Input, integer NTAB, the number of data points XTAB.
%
%    Input, real XTAB(NTAB), the X values upon which the
%    Lagrange basis polynomial is to be based.
%
%    Output, real DIFTAB(NTAB), the divided difference table
%    for the IVAL-th Lagrange basis polynomial.
%

%
%  Check IVAL.
%
  if ( ival < 1 | ntab < ival )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DIF_BASIS_I - Fatal error!\n' );
    fprintf ( 1, '  IVAL must be between 1 and %d\n', ntab );
    fprintf ( 1, '  but your value is %d\n', ival );
    error ( 'DIF_BASIS_I - Fatal error!' );
  end
%
%  Initialize DIFTAB to Delta(I,J).
%
  diftab(1:ntab) = 0.0;
  diftab(ival) = 1.0;
%
%  Compute the IVAL-th Lagrange basis polynomial.
%
  diftab = data_to_dif ( ntab, xtab, diftab );

  return
end

author: johannesgerer


text:
function [B_detail, B_smooth] = train_choose_detail(training_path, unit, col)
% choosing image patches which contain more detail

fileFolder=fullfile(training_path);
dirOutput=dir(fullfile(fileFolder,'*'));
num = length(dirOutput);
% fileNames={dirOutput.name}'; % all names
B = [];
for i = 3:num
    img = imread([training_path,dirOutput(i).name]);
    if size(img,3)>1
        img = rgb2gray(img);
    end

    [t1,t2] = size(img);
    img = img(1:floor(t1/unit)*unit, 1:floor(t2/unit)*unit);

    img = im2double(img);
    b = im2col(img, [unit, unit], 'distinct');
    B = cat(2, B, b);
end
% calculating SD valuee of each column B
B_mean = mean(B,1);
B_mean = repmat(B_mean, unit*unit, 1);
B_sd = sqrt(sum((B - B_mean).*(B - B_mean),1));
% choose >0.5
B_sd(find(B_sd<=0.5)) = 0;
B_sd(find(B_sd>0.5)) = 1;
B_mask = repmat(B_sd, unit*unit, 1);
B_smooth = B.*(1-B_mask);
B_detail = B.*B_mask;
B_smooth(:,all(B_smooth==0,1))=[];
B_detail(:,all(B_detail==0,1))=[];

[m1,n1] = size(B_detail);
[m2,n2] = size(B_smooth);

end

author: thfylsty


text:
function pot = convert_to_pot(CPD, pot_type, domain, evidence)
% CONVERT_TO_POT Convert a gmux CPD to a Gaussian potential
% pot = convert_to_pot(CPD, pot_type, domain, evidence)

switch pot_type
 case {'d', 'u', 'cg', 'scg'},
  error(['can''t convert gmux to potential of type ' pot_type])

 case {'c','g'},
  % We create a large weight matrix with zeros in all blocks corresponding
  % to the non-chosen parents, since they are effectively disconnected.
  % The chosen parent is determined by the value, m,  of the discrete parent.
  % Thus the potential is as large as the whole family.
  ps = domain(1:end-1);
  dps = ps(CPD.dps); % CPD.dps is an index, not a node number (because of param tying)
  cps = ps(CPD.cps);
  m = evidence{dps};
  if isempty(m)
    error('gmux node must have observed discrete parent')
  end
  bs = CPD.sizes(CPD.cps);
  b = block(m, bs);
  sum_cpsz = sum(CPD.sizes(CPD.cps));
  selfsz = CPD.sizes(end);
  W = zeros(selfsz, sum_cpsz);
  W(:,b) = CPD.weights(:,:,m);

  ns = zeros(1, max(domain));
  ns(domain) = CPD.sizes;
  self = domain(end);
  cdom = [cps(:)' self];
  pot = linear_gaussian_to_cpot(CPD.mean(:,m), CPD.cov(:,:,m), W, domain, ns, cdom, evidence);

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

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function self = correct(self)
     ref = 0.707; % note that median, rms and std all give same value on x=sin(0:pi/1000:2*pi)
     for c=1:numel(self)
        if strcmp(self(c).measure, 'max')
            self(c).data = self(c).data * ref;
        end
        if strcmp(self(c).measure, '68')
            self(c).data = self(c).data/0.8761 * ref;
        end
        if strcmp(self(c).measure, 'mean')
            self(c).data = self(c).data/0.6363 * ref;
        end
     end
end

author: geoscience-community-codes

repo: GISMO


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

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

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

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

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

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

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

% for Channel D
chD_wdth   = [   10     10     10     10     10     10     10     10     10     20     30     30     30     40     50     50     50     50];

% define imag. component
j = sqrt(-1);

% Initialize channel model
if (Clk==0)
    % Initialize PER test
    snr_idx = 1;
    % Load Channel D coeff's
    load tgn_d_nlos;
    % Setup channel coeff's for PER test
    H_cohere_prd = ceil(size(H, 4)/100);  % ch. mat samples in one coherence time (100 simulated)
    H_idx_inc = round(H_cohere_prd / numSampPerCohere);
    % Setup indices for channel coeff's
    num_taps = size(H, 3);
    H_idx_ini = 100; %100; %2000;  % initial index for channel matrices
    H_idx = H_idx_ini;
    H_cnt = 1;

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

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

 % Reuse channel coeff's for each SNR tested...
 if (mod(Clk,numPkts)==(numPkts-1))  H_idx = H_idx_ini;  end
end

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

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

    % Apply channel filter
    H_dline_idx = 1;
    H_out = zeros(Nrx,1);
    for m=1:num_taps
if 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Apply tap to inputs (considering tap width)
        H_samp = zeros(Ntx, 1);
        for k = 1:(chD_wdth(m)/10)
            % 10 ns (100 MHz) sampling
            H_samp = H_samp + H_dline(:,H_dline_idx);
            H_dline_idx = H_dline_idx + 1;
        end
        H_out = H_out + H_tap(:,:,m)*H_samp;
else %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Apply tap to inputs (considering tap width)
        for k = 1:(chD_wdth(m)/10)
  if 1 % 10 ns (100 MHz) sampling
            H_out = H_out + H_tap(:,:,m)*H_dline(:,H_dline_idx);
  else % 50 ns (20 MHz) sampling
            if mod(H_dline_idx,5)==1
                H_out = H_out + H_tap(:,:,m)*H_dline(:,1+floor(H_dline_idx/5));
            end
  end
            H_dline_idx = H_dline_idx + 1;
        end
end  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end

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

% Add AWGN to Rx streams (reduce SNR by 5 dB, due to upsampling)
snr_val = PER_snr(snr_idx);
for i=1:Nrx
    ch_out(i,:) = awgn(ch_out(i,:), snr_val-5, 'measured');
end

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

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

author: Sable

repo: mcbench-benchmarks


text:
function value = nco_abscissa ( order, i )

%*****************************************************************************80
%
%% NCO_ABSCISSA returns the I-th abscissa for the Newton Cotes open rule.
%
%  Discussion:
%
%    Our convention is that the abscissas are numbered from left to
%    right.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 May 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ORDER, the order of the rule.
%    1 <= ORDER.
%
%    Input, integer I, the index of the desired abscissa.
%    1 <= I <= ORDER.
%
%    Output, real VALUE, the value of the I-th
%    abscissa in the Newton Cotes open rule of order ORDER.
%
  x_min = -1.0;
  x_max = +1.0;

  if ( order < 1 )
    value = - Inf;
    return
  end

  if ( i < 1 | order < i )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'NCO_ABSCISSA - Fatal error!\n' );
    fprintf ( 1, '  1 <= I <= ORDER is required.\n' );
    error ( 'NCO_ABSCISSA - Fatal error!' );
  end

  value = ( ( order - i + 1 ) * x_min   ...
          + (         i     ) * x_max ) ...
          / ( order     + 1 );

  return
end

repo: jburkardt-m


text:
function Cpp_Params = Flatten_Cpp_Params(P_Cpp,Data)

%FUNCTION:
%   This function takes a state struct S, and collapses it to a column
%   vector s. It relies on information about the struct format that is
%   stored in P.StateStruct, which is created by the
%   Define_StateStruct_Data.m function.
%
%INPUTS:
%   P = the struct the is being used ito pass parameters
%   Data = contains information about P
%
%OUTPUTS:
%   param_vec = a column vector with all of the data in P
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Names = Data.Names;
Sizes = Data.Sizes;
Idx = Data.Idx;

Length = Idx(end,2);   %How many elements in the entire state
Cpp_Params = zeros(Length,1);   %Total size of the state struct
N = length(Names);   %total number of fields;
L = Sizes(:,1).*Sizes(:,2);   %Length of each field;

for i=1:N
   Cpp_Params(Idx(i,1):Idx(i,2)) = reshape(P_Cpp.(Names{i})',L(i),1);  %Transpose accounts for C++ array definition difference
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function C = cross_product(A,B)

% Calculates the cross product C of the 3-vectors A and B

C = [A(2)*B(3)-A(3)*B(2);  A(3)*B(1)-A(1)*B(3);  A(1)*B(2)-A(2)*B(1)];

author: InverseTampere

repo: TreeQSM


text:
%This directory contains the files that are required to calibrate a
%BBCI processing model, to save it, and to apply it to continuously
%incoming data.
%
%BBCI_CALIBRATE - Establish a classifier based on calibration data
%BBCI_CALIBRATE_STRUCTURES - Help file describing the data structures
%
%BBCI_SAVE - Save BBCI classifier (and all information required for
%            online operation in a file.
%
%BBCI_APPLY - Apply BBCI classifier to continuously acquired data
%BBCI_APPLY_STRUCTRES - Help file describing the data structures

author: bbci


text:
%stepx   outlines step domain
%   IFISS scriptfile: DJS; 28 February 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
if exist('x','var')==0,
   error('You need to set up a step-shaped domain first!'),
end
varsteplen = max(x);
hold on
plot([-1,0],[0,0],'-k');
plot([0,0],[0,-1],'-k');
plot([0,varsteplen],[-1,-1],'-k');
plot([varsteplen,varsteplen],[-1,1],'-k');
plot([varsteplen,-1],[1,1],'-k');
plot([-1,-1],[1,0],'-k');
hold off


text:
function [res,R,p,sigma,Ms] = TransformK(p,rho,sigma,Ms)
%TRANSFORMK   Transformation res+R+sum(ps) = rho+sum(p)
%
%   [res,R,ps] = TransformK(p,rho)
%
%Input vector p may be single or double precision.
%
%Implements Algorithm 6.2 from
%  S.M. Rump, T. Ogita, S. Oishi: Accurate Floating-point Summation II:
%    Sign, K-fold Faithful and Rounding to Nearest, Siam J. Sci. Comput.,
%    31(2):1269-1302, 2008.
%Requires (4m+3)n flops for m executions of repeat-until loop if
%  sigma,Ms not given, otherwise (4m+1)n flops.
%
%Reference implementation! Slow due to interpretation!
%

% written  03/03/07     S.M. Rump
%

  if nargin==4                          % sigma and Ms already known
    kPhi = 0;
    [tau1,tau2,p,sigma,Ms] = Transform(p,rho,kPhi,sigma,Ms);
  else                                  % standard call
    [tau1,tau2,p,sigma,Ms] = Transform(p,rho);
  end
  res = tau1 + ( tau2 + sum(p) );
  R = tau2 - ( res - tau1 );

repo: OpenMAS


text:
classdef TuneAGC < adi.common.DebugAttribute & adi.common.RegisterReadWrite
    properties (Nontunable, Hidden)
        CustomAGC = 0;

        AttackDelay = 1;
        PeakOverloadWaitTime = 10;
        AGCLockLevel = 10;
        DecStepSizeFullTableCase3 = 3;
        ADCLargeOverloadThresh = 58;
        ADCSmallOverloadThresh = 47;
        DecStepSizeFullTableCase2 = 3;
        DecStepSizeFullTableCase1 = 3;
        LargeLMTOverloadThresh = 35;
        SmallLMTOverloadThresh = 25;
        SettlingDelay = 3;
        EnergyLostThresh = 3;
        LowPowerThresh = 15;
        IncrementGainStep
        FAGCLockLevelGainIncreaseUpperLimit  = 7;
        FAGCLPThreshIncrementTime = 3;
        DecPowMeasurementDuration = 16;
    end

    properties (Constant, Hidden, Access = private)
        % Register addresses in hexadecimal
        AttackDelay_Reg = '022';
        PeakOverloadWaitTime_Reg = '0FE';
        AGCLockLevel_Reg = '101';
        DecStepSizeFullTableCase3_Reg = '103';
        ADCSmallOverloadThresh_Reg = '104';
        ADCLargeOverloadThresh_Reg = '105';
        DecStepSizeFullTableCase2_Reg = '106';
        DecStepSizeFullTableCase1_Reg = '106';
        LargeLMTOverloadThresh_Reg = '108';
        SmallLMTOverloadThresh_Reg = '107';
        SettlingDelay_Reg = '111';
        EnergyLostThresh_Reg = '112';
        LowPowerThresh_Reg = '114';
        IncrementGainStep_Reg = '117';
        FAGCLockLevelGainIncreaseUpperLimit_Reg = '118';
        FAGCLPThreshIncrementTime_Reg = '11B';
        DecPowMeasurementDuration_Reg = '15C';

        % Register mask in binary
        AttackDelay_Mask = '11000000';
        PeakOverloadWaitTime_Mask = '11100000';
        AGCLockLevel_Mask = '10000000';
        DecStepSizeFullTableCase3_Mask = '11100011';
        DecStepSizeFullTableCase2_Mask = '10001111';
        DecStepSizeFullTableCase1_Mask = '11110000';
        LargeLMTOverloadThresh_Mask = '11000000';
        SmallLMTOverloadThresh_Mask = '11000000';
        SettlingDelay_Mask = '11100000';
        EnergyLostThresh_Mask = '11000000';
        LowPowerThresh_Mask = '10000000';
        IncrementGainStep_Mask = '00011111';
        FAGCLockLevelGainIncreaseUpperLimit_Mask = '11000000';
        DecPowMeasurementDuration_Mask = '11110000';

        % Bit-shifts to be applied
        DecStepSizeFullTableCase3_BitShift = 2;
        DecStepSizeFullTableCase2_BitShift = 4;
        IncrementGainStep_BitShift = 5;
    end

    methods
        function set.AttackDelay(obj, value)
            validateattributes( value, { 'double','single', 'uint32' }, ...
                { 'real', 'nonnegative','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',0,'<=',63}, ...
                '', 'AttackDelay');
            obj.AttackDelay = value;
            if obj.ConnectedToDevice
                obj.setRegister(value, obj.AttackDelay_Reg, obj.AttackDelay_Mask);
            end
        end
        function set.PeakOverloadWaitTime(obj, value)
            validateattributes( value, { 'double','single', 'uint32' }, ...
                { 'real', 'nonnegative','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',0,'<=',31}, ...
                '', 'PeakOverloadWaitTime');
            obj.PeakOverloadWaitTime = value;
            if obj.ConnectedToDevice
                obj.setRegister(value, obj.PeakOverloadWaitTime_Reg, obj.PeakOverloadWaitTime_Mask);
            end
        end
        function set.AGCLockLevel(obj, value)
            validateattributes( value, { 'double','single', 'uint32' }, ...
                { 'real', 'nonnegative','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',0,'<=',127}, ...
                '', 'AGCLockLevel');
            obj.AGCLockLevel = value;
            if obj.ConnectedToDevice
                obj.setRegister(value, obj.AGCLockLevel_Reg, obj.AGCLockLevel_Mask);
            end
        end
        function set.DecStepSizeFullTableCase3(obj, value)
            validateattributes( value, { 'double','single', 'uint32' }, ...
                { 'real', 'nonnegative','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',0,'<=',7}, ...
                '', 'DecStepSizeFullTableCase3');
            obj.DecStepSizeFullTableCase3 = value;
            if obj.ConnectedToDevice
                obj.setRegister(value, obj.DecStepSizeFullTableCase3_Reg, obj.DecStepSizeFullTableCase3_Mask, obj.DecStepSizeFullTableCase3_BitShift);
            end
        end
        function set.ADCLargeOverloadThresh(obj, value)
            validateattributes( value, { 'double','single', 'uint32' }, ...
                { 'real', 'nonnegative','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',0,'<=',255}, ...
                '', 'ADCLargeOverloadThresh');
            obj.ADCLargeOverloadThresh = value;
            if obj.ConnectedToDevice
                obj.setDebugAttributeLongLong('adi,gc-adc-large-overload-thresh',value);
            end
        end
        function set.ADCSmallOverloadThresh(obj, value)
            validateattributes( value, { 'double','single', 'uint32' }, ...
                { 'real', 'nonnegative','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',0,'<=',obj.ADCLargeOverloadThresh}, ...
                '', 'ADCSmallOverloadThresh');
            obj.ADCSmallOverloadThresh = value;
            if obj.ConnectedToDevice
                obj.setDebugAttributeLongLong('adi,gc-adc-small-overload-thresh',value);
            end
        end
        function set.DecStepSizeFullTableCase2(obj, value)
            validateattributes( value, { 'double','single', 'uint32' }, ...
                { 'real', 'nonnegative','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',0,'<=',7}, ...
                '', 'DecStepSizeFullTableCase2');
            obj.DecStepSizeFullTableCase2 = value;
            if obj.ConnectedToDevice
                obj.setRegister(value, obj.DecStepSizeFullTableCase2_Reg, obj.DecStepSizeFullTableCase2_Mask, obj.DecStepSizeFullTableCase2_BitShift);
            end
        end
        function set.DecStepSizeFullTableCase1(obj, value)
            validateattributes( value, { 'double','single', 'uint32' }, ...
                { 'real', 'nonnegative','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',0,'<=',15}, ...
                '', 'DecStepSizeFullTableCase1');
            obj.DecStepSizeFullTableCase1 = value;
            if obj.ConnectedToDevice
                obj.setRegister(value, obj.DecStepSizeFullTableCase1_Reg, obj.DecStepSizeFullTableCase1_Mask);
            end
        end
        function set.LargeLMTOverloadThresh(obj, value)
            validateattributes( value, { 'double','single', 'uint32' }, ...
                { 'real', 'nonnegative','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',0,'<=',63}, ...
                '', 'LargeLMTOverloadThresh');
            obj.LargeLMTOverloadThresh = value;
            if obj.ConnectedToDevice
                obj.setRegister(value, obj.LargeLMTOverloadThresh_Reg, obj.LargeLMTOverloadThresh_Mask);
            end
        end
        function set.SmallLMTOverloadThresh(obj, value)
            validateattributes( value, { 'double','single', 'uint32' }, ...
                { 'real', 'nonnegative','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',0,'<=',obj.LargeLMTOverloadThresh}, ...
                '', 'SmallLMTOverloadThresh');
            obj.SmallLMTOverloadThresh = value;
            if obj.ConnectedToDevice
                obj.setRegister(value, obj.SmallLMTOverloadThresh_Reg, obj.SmallLMTOverloadThresh_Mask);
            end
        end
        function set.SettlingDelay(obj, value)
            validateattributes( value, { 'double','single', 'uint32' }, ...
                { 'real', 'nonnegative','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',0,'<=',31}, ...
                '', 'SettlingDelay');
            obj.SettlingDelay = value;
            if obj.ConnectedToDevice
                obj.setRegister(value, obj.SettlingDelay_Reg, obj.SettlingDelay_Mask);
            end
        end
        function set.EnergyLostThresh(obj, value)
            validateattributes( value, { 'double','single', 'uint32' }, ...
                { 'real', 'nonnegative','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',0,'<=',63}, ...
                '', 'SettlingDelay');
            obj.EnergyLostThresh = value;
            if obj.ConnectedToDevice
                obj.setRegister(value, obj.EnergyLostThresh_Reg, obj.EnergyLostThresh_Mask);
            end
        end
        function set.LowPowerThresh(obj, value)
            validateattributes( value, { 'double','single', 'uint32' }, ...
                { 'real', 'nonnegative','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',0,'<=',63}, ...
                '', 'LowPowerThresh');
            obj.LowPowerThresh = value;
            if obj.ConnectedToDevice
                obj.setDebugAttributeLongLong('adi,gc-low-power-thresh',value);
            end
        end
        function set.IncrementGainStep(obj, value)
            validateattributes( value, { 'double','single', 'uint32' }, ...
                { 'real', 'nonnegative','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',0,'<=',7}, ...
                '', 'IncrementGainStep');
            obj.IncrementGainStep = value;
            if obj.ConnectedToDevice
                obj.setRegister(value, obj.IncrementGainStep_Reg, obj.IncrementGainStep_Mask, obj.IncrementGainStep_BitShift);
            end
        end
        function set.FAGCLockLevelGainIncreaseUpperLimit(obj, value)
            validateattributes( value, { 'double','single', 'uint32' }, ...
                { 'real', 'nonnegative','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',0,'<=',63}, ...
                '', 'FAGCLockLevelGainIncreaseUpperLimit');
            obj.FAGCLockLevelGainIncreaseUpperLimit = value;
            if obj.ConnectedToDevice
                obj.setDebugAttributeLongLong('adi,fagc-lock-level-gain-increase-upper-limit',value);
            end
        end
        function set.FAGCLPThreshIncrementTime(obj, value)
            validateattributes( value, { 'double','single', 'uint32' }, ...
                { 'real', 'nonnegative','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',0,'<=',255}, ...
                '', 'FAGCLPThreshIncrementTime');
            obj.FAGCLPThreshIncrementTime = value;
            if obj.ConnectedToDevice
                obj.setDebugAttributeLongLong('adi,fagc-lp-thresh-increment-time',value);
            end
        end
        function set.DecPowMeasurementDuration(obj, value)
            validateattributes( value, { 'double','single', 'uint32' }, ...
                { 'real', 'nonnegative','scalar', 'finite', 'nonnan', 'nonempty','integer','>=',0,'<=',15}, ...
                '', 'DecPowMeasurementDuration');
            obj.DecPowMeasurementDuration = value;
            if obj.ConnectedToDevice
                obj.setRegister(value, obj.DecPowMeasurementDuration_Reg, obj.DecPowMeasurementDuration_Mask);
            end
        end
        function WriteDebugAttributes(obj)
            if obj.ConnectedToDevice
                obj.setDebugAttributeLongLong('adi,gc-adc-large-overload-thresh',obj.ADCLargeOverloadThresh);
                obj.setDebugAttributeLongLong('adi,gc-adc-small-overload-thresh',obj.ADCSmallOverloadThresh);
                obj.setDebugAttributeLongLong('adi,gc-low-power-thresh',obj.LowPowerThresh);
                obj.setDebugAttributeLongLong('adi,fagc-lock-level-gain-increase-upper-limit',obj.FAGCLockLevelGainIncreaseUpperLimit);
                obj.setDebugAttributeLongLong('adi,fagc-lp-thresh-increment-time',obj.FAGCLPThreshIncrementTime);
            end
        end
        function WriteToRegisters(obj)
            if obj.ConnectedToDevice
                obj.setRegister(obj.AttackDelay, obj.AttackDelay_Reg, obj.AttackDelay_Mask);
                obj.setRegister(obj.PeakOverloadWaitTime, obj.PeakOverloadWaitTime_Reg, obj.PeakOverloadWaitTime_Mask);
                obj.setRegister(obj.AGCLockLevel, obj.AGCLockLevel_Reg, obj.AGCLockLevel_Mask);
                obj.setRegister(obj.DecStepSizeFullTableCase3, obj.DecStepSizeFullTableCase3_Reg, obj.DecStepSizeFullTableCase3_Mask, obj.DecStepSizeFullTableCase3_BitShift);
                obj.setRegister(obj.DecStepSizeFullTableCase2, obj.DecStepSizeFullTableCase2_Reg, obj.DecStepSizeFullTableCase2_Mask, obj.DecStepSizeFullTableCase2_BitShift);
                obj.setRegister(obj.DecStepSizeFullTableCase1, obj.DecStepSizeFullTableCase1_Reg, obj.DecStepSizeFullTableCase1_Mask);
                obj.setRegister(obj.LargeLMTOverloadThresh, obj.LargeLMTOverloadThresh_Reg, obj.LargeLMTOverloadThresh_Mask);
                obj.setRegister(obj.SmallLMTOverloadThresh, obj.SmallLMTOverloadThresh_Reg, obj.SmallLMTOverloadThresh_Mask);
                obj.setRegister(obj.SettlingDelay, obj.SettlingDelay_Reg, obj.SettlingDelay_Mask);
                obj.setRegister(obj.EnergyLostThresh, obj.EnergyLostThresh_Reg, obj.EnergyLostThresh_Mask);
                obj.setRegister(obj.IncrementGainStep, obj.IncrementGainStep_Reg, obj.IncrementGainStep_Mask, obj.IncrementGainStep_BitShift);
                obj.setRegister(obj.DecPowMeasurementDuration, obj.DecPowMeasurementDuration_Reg, obj.DecPowMeasurementDuration_Mask);
            end
        end
        function value = ReadFromRegister(obj, prop_name)
            if obj.ConnectedToDevice
                switch prop_name
                    case 'AttackDelay'
                        value = obj.getRegister(obj.AttackDelay_Reg, obj.AttackDelay_Mask);
                    case 'PeakOverloadWaitTime'
                        value = obj.getRegister(obj.PeakOverloadWaitTime_Reg, obj.PeakOverloadWaitTime_Mask);
                    case 'AGCLockLevel'
                        value = obj.getRegister(obj.AGCLockLevel_Reg, obj.AGCLockLevel_Mask);
                    case 'DecStepSizeFullTableCase3'
                        value = obj.getRegister(obj.DecStepSizeFullTableCase3_Reg, obj.DecStepSizeFullTableCase3_Mask, obj.DecStepSizeFullTableCase3_BitShift);
                    case 'ADCSmallOverloadThresh'
                        value = obj.getRegister(obj.ADCSmallOverloadThresh_Reg);
                    case 'ADCLargeOverloadThresh'
                        value = obj.getRegister(obj.ADCLargeOverloadThresh_Reg);
                    case 'DecStepSizeFullTableCase2'
                        value = obj.getRegister(obj.DecStepSizeFullTableCase2_Reg, obj.DecStepSizeFullTableCase2_Mask, obj.DecStepSizeFullTableCase2_BitShift);
                    case 'DecStepSizeFullTableCase1'
                        value = obj.getRegister(obj.DecStepSizeFullTableCase1_Reg, obj.DecStepSizeFullTableCase1_Mask);
                    case 'LargeLMTOverloadThresh'
                        value = obj.getRegister(obj.LargeLMTOverloadThresh_Reg, obj.LargeLMTOverloadThresh_Mask);
                    case 'SmallLMTOverloadThresh'
                        value = obj.getRegister(obj.SmallLMTOverloadThresh_Reg, obj.SmallLMTOverloadThresh_Mask);
                    case 'SettlingDelay'
                        value = obj.getRegister(obj.SettlingDelay_Reg, obj.SettlingDelay_Mask);
                    case 'EnergyLostThresh'
                        value = obj.getRegister(obj.EnergyLostThresh_Reg, obj.EnergyLostThresh_Mask);
                    case 'LowPowerThresh'
                        value = obj.getRegister(obj.LowPowerThresh_Reg, obj.LowPowerThresh_Mask);
                    case 'IncrementGainStep'
                        value = obj.getRegister(obj.IncrementGainStep_Reg, obj.IncrementGainStep_Mask, obj.IncrementGainStep_BitShift);
                    case 'FAGCLockLevelGainIncreaseUpperLimit'
                        value = obj.getRegister(obj.FAGCLockLevelGainIncreaseUpperLimit_Reg, obj.FAGCLockLevelGainIncreaseUpperLimit_Mask);
                    case 'FAGCLPThreshIncrementTime'
                        value = obj.getRegister(obj.FAGCLPThreshIncrementTime_Reg);
                    case 'DecPowMeasurementDuration'
                        value = obj.getRegister(obj.DecPowMeasurementDuration_Reg, obj.DecPowMeasurementDuration_Mask);
                    otherwise
                        error('Attempted to read unknown property %s\n', prop_name);
                end
            end
        end
    end
end

author: analogdevicesinc


text:
function example
% Two-state model of gene expression
%
% Reaction network:
%   0 -> mRNA
%   mRNA -> mRNA + protein
%   mRNA -> 0
%   protein -> 0

tspan = [0, 10000]; %seconds
x0 = [0, 0]; %mRNA, protein
stoich_matrix = [ 1  0  ; %transcription
                  0  1  ; %translation
                 -1  0  ; %mRNA degradation
                  0 -1 ]; %protein degradation

% Rate constants
p.kR = 0.1;%0.01;
p.kP = 0.1;%1;
p.gR = 0.1;
p.gP = 0.002;

% Run simulation
%[t,x] = directMethod(stoich_matrix, @propensities_2state, tspan, x0, p);
[t,x] = firstReactionMethod(stoich_matrix, @propensities_2state, tspan, x0, p);

% Plot time course
figure(gcf);
stairs(t,x);
set(gca,'XLim',tspan);
xlabel('time (s)');
ylabel('molecules');
legend({'mRNA','protein'});

end

function a = propensities_2state(x, p)
mRNA    = x(1);
protein = x(2);

a = [p.kR;
     p.kP*mRNA;
	 p.gR*mRNA;
     p.gP*protein];
end


text:
function [W,E,U] = skeleton_extraction(V,F)
  % SKELETON_EXTRACTION Compute the "skeleton" of a surface mesh following "Skeleton
  % Extraction by Mesh Contraction" [Au et al. 2008]. The final combinatorially
  % simplification is a bit different than the QSlim-based approach in [Au et
  % al. 2008].
  %
  % [W,E,U] = skeletonization(V,F)
  %
  % Inputs:
  %   V  #V by 3 list of mesh vertex positions
  %   F  #F by 3 list of mesh indices
  % Outputs:
  %   W  #W by 3 list of skeleton vertex positions
  %   E  #E by 2 list of skeleton edges into W
  %   U  #V by 3 list of skeleton vertex positions, before thinning
  %

  delta = 1e-3*max(max(V)-min(V))^2;
  U = V;

  viz = false;
  if viz
      tsurf(F,V,'FaceAlpha',0.1,'EdgeColor','k','FaceColor',[0.5 0.5 0.5],'EdgeAlpha',0.5);
      hold on;
      t = tsurf( ...
          bsxfun(@plus,size(F,1)*(0:2),(1:size(F,1))'), ...
          U(F,:),'FaceColor',[0.7 0.8 1.0],'EdgeColor','k','EdgeAlpha',0.5, ...
          'SpecularStrength',0.1,'FaceLighting','phong');
      tr = tsurf([1 1 1],U,'EdgeColor','r','LineWidth',2);
      hold off;
      axis equal;
      view(-62,10);
      camproj('persp');
      light('Position',[-100.0,1.0,-1.0],'Style','infinite');
      set(gca,'Visible','off');
      set(gcf,'Color','w');
      drawnow;
  end

  b = [];

  iter = 1;
  h = avgedge(V,F);

  while true
    L = cotmatrix(U,F);
    M = massmatrix(U,F,'barycentric');
    b = union(find(any(abs(L)>1e+5,2)),b);
    U_prev = U;
    U = min_quad_with_fixed(-(M-delta*L),2*M*U,b,U(b,:));
    %U = U/sqrt(sum(doublearea(U,F))*0.5);
    %U = bsxfun(@minus,U,centroid(U,F));
    %D = max(internalangles(U,F),[],2)>0.9*pi & ...
    %  doublearea(U,F)<1e-5;
    D = all(ismember(F,b),2);
    if viz
        set(t,'Vertices',U(F,:));
        set(tr,'Faces',F(D,:),'Vertices',U);
        xlabel(sprintf('%d',numel(b)));
        drawnow;
    end
    if max(max(abs(U-U_prev))) < 1e-2*h
      break;
    end
    iter =iter + 1;
  end

  W = U;
  A = adjacency_matrix(F) + speye(size(V,1));
  EC = adjacency_edge_cost_matrix(U,F);
  phi = zeros(size(W,1),1);
  marked = false(size(W,1),1);
  E = [];
  % This is very slow for large meshes...
  for wi = 1:size(W,1)
    % this is a lot faster than find(A(ui,:))
    if marked(wi)
      continue;
    end
    pI = A(:,wi)>0;
    I = false(size(V,1),1);
    N = [];
    while true
      pN = find(pI);
      if numel(pN) <= 1
        break;
      end
      Wi = W(pN,:);
      [pC,pS,L] = pca(Wi);
      phi(wi) = L(1)/sum(L);
      if phi(wi) < 0.99
        break;
      end
      I = pI;
      N = pN;
      C = pC;
      S = pS;

      if viz
          tsurf(F,W);
          hold on;
          scatter3(U(I,1),U(I,2),U(I,3),'SizeData',50);
          hold off;
          title(sprintf('%g',phi(wi)));
          axis equal;
          drawnow;
      end

      %pI = (A*I)>0;
      CCI = sum(A(I,:))';
      ECI = minnz(EC(I,:))';
      ECI = ECI./CCI;
      ECI(I) = 0;
      [~,j] = minnz(ECI);
      pI = I;
      pI(j) = 1;

      if isequal(pI,I)
        break;
      end

    end

    if ~isempty(N)
      [~,O] = sort(S(:,1));
      O = N(O);
      A(I,:) = 0;
      A(:,I) = 0;
      %tsurf(F,W);
      %hold on;
      %scatter3(U(I,1),U(I,2),U(I,3),'SizeData',25);
      %scatter3(U(marked,1),U(marked,2),U(marked,3),'SizeData',50);
      %plot_edges(U,E,'b','LineWidth',3);
      %plot_edges(U,[O(1) O(end)],'r','LineWidth',3);
      %hold off;
      %axis equal;
      %drawnow
      marked(N) = true;
      E = [E;O(1) O(end)];
    end

    %if phi<0.4
    %  continue;
    %end
    %Wi = bsxfun(@plus,bsxfun(@times,S(:,1),C(:,1)'),mean(Wi));
    %d = sqrt(sum(bsxfun(@minus,Wi,Wi(1,:)).^2,2));
    %H = max(d);
    %d = d/H;
    %f = 2.*d.^3 - 3.*d.^2 + 1;
    %f = f/sum(f);
    %W(wi,:) = sum(bsxfun(@times,f,Wi));
  end
  %plot_edges(W,E,'LineWidth',2)

  count = ones(size(W,1),1);
  while true
    [W,IM] = remove_unreferenced(W,E);
    count(IM) = count;
    count = count(1:size(W,1));
    E = IM(E);
    %plot_edges(W,E,'LineWidth',2)
    %axis equal;
    %drawnow;
    if size(E,1) == 0
      break;
    end
    A = adjacency_matrix(E)+speye(size(W,1));
    [ncc,C] = conncomp(A);
    if ncc == 1
      break;
    end
    D = pdist2(W,W);
    D(A>0) = inf;
    [D,I] = min(D);
    [~,j] = min(D);
    i = I(j);
    % running average
    W(i,:) = (count(i)*W(i,:)+W(j,:))/(count(i)+1);
    count(i) = count(i) + 1;
    E(E==j) = i;
    E = unique(sort(E,2),'rows');
  end

  %G = F;
  %E = edges(G);
  %
  %A = U(E(:,2),:)-U(E(:,1),:);
  %B = cross(A,U(E(:,1),:),2);
  %UH = [U ones(size(U,1),1)];
  %
  %ne = size(E,1);
  %n = size(U,1);
  %%K = sparse( ...
  %%  repmat(reshape(repmat(1:ne*3,1,3),ne,9),2,1), ...
  %%  bsxfun(@plus,E(:),n*[1 2 0 2 0 1 3 3 3]), ...
  %%  repmat([-A(:,[3 1 2]) A(:,[2 3 1]) -B],2,1), ...
  %%  ne*3,n*4);
  %%Q = K'*K;
  %
  %%K = zeros(3,4,size(E,1));
  %Q = zeros(4,4,size(U,1));
  %%Q = sparse(n*4,n*4);
  %% This can at least be sped up using http://www.alecjacobson.com/weblog/?p=4186
  %for ei = 1:size(E,1)
  %  K = [ ...
  %    0        -A(ei,3)  A(ei,2) -B(ei,1); ...
  %    A(ei,3)  0        -A(ei,1) -B(ei,2); ...
  %    -A(ei,2) A(ei,1)  0        -B(ei,3)];
  %  K2 = K'*K;
  %  Q(:,:,E(ei,1)) = Q(:,:,E(ei,1)) + K2;
  %  Q(:,:,E(ei,2)) = Q(:,:,E(ei,2)) + K2;
  %end
  %% Cost per vertex
  %C = zeros(size(U,1),1);
  %for vi = 1:size(U,1)
  %  C(vi) = UH(vi,:) * (Q(:,:,vi) * (UH(vi,:)'));
  %end
  %
  %
  %A = sparse(repmat((1:ne)',1,2),E,1,ne,n);
  %EL = sqrt(sum(U(E(:,2),:)-U(E(:,1),:),2));
  %% This is a terrible O(n²) implemenation. The energy also doesn't match the [Au
  %% et al. 2008] paper. I don't understand their description or necessarily agree
  %% with the motivation.
  %J = 1:size(U,1);
  %while true
  %  % Edge-vertex incidence matrix
  %  ne = size(E,1);
  %  % Edge lengths summed at vertices
  %  EC = repmat(A*C,2,1) + 0.001*repmat(EL,2,1);
  %  % pick an edge
  %  [~,ei] = min(EC);
  %  e = E(mod(ei-1,ne)+1,:);
  %  if ei>ne
  %    e = fliplr(e);
  %  end
  %  % collapse edge
  %  J(e(2)) = e(1);
  %  J = J(J);
  %  Q(:,:,e(1)) = Q(:,:,e(1)) + Q(:,:,e(2));
  %  C(e(1)) = UH(e(1),:) * (Q(:,:,e(1)) * (UH(e(1),:)'));
  %  G = J(G);
  %  E = J(E);
  %
  %  A(:,e(1)) = A(:,e(1)) + A(:,e(2));
  %  A(:,e(2)) = 0;
  %  ndeg = E(:,1)~=E(:,2);
  %  E = E(ndeg,:);
  %  EL = EL(ndeg,:);
  %  A = A(ndeg,:);
  %  G = G(G(:,1)~=G(:,2)&G(:,2)~=G(:,3)&G(:,3)~=G(:,1),:);
  %  %G = unique(sort(G,2),'rows');
  %  if isempty(G)
  %    break;
  %  end
  %
  %  change = any(E==e(1),2);
  %  Echange = unique(sort(E(change,:),2),'rows');
  %  E = [E(~change,:);Echange];
  %  nc = size(Echange,1);
  %  A = [ ...
  %    A(~change,:); ...
  %    sparse(repmat((1:nc)',1,2),Echange,1,nc,n)];
  %  EL = [EL(~change);sqrt(sum(U(Echange(:,2),:)-U(Echange(:,1),:),2))];
  %
  %  if size(E,1)<50
  %    tsurf(G,U,'EdgeColor','r');
  %    hold on;
  %    tsurf([E E(:,1)],U,'EdgeColor','k');
  %    hold off;
  %    title(sprintf('%d,%d',size(G,1),size(E,1)));
  %    axis equal;
  %    drawnow
  %  end
  %end
  %M = sparse(J,1:n,diag(massmatrix(V,F)),n,n);
  %M(J,:) = M(J,:)*diag(1./sum(M(J,:),2));
  %U = M*V;
  %[U,IM] = remove_unreferenced(U,E);
  %E = IM(E);
  %plot_edges(U,E,'LineWidth',2);
  %hold on;
  %tsurf(F,V,'CData',J,'FaceAlpha',0.2,'EdgeAlpha',0.2);
  %hold off;
end

author: alecjacobson

repo: gptoolbox


text:
classdef (Abstract) MimStorageClass  < handle

    methods (Abstract)
        [metaData, dataStream] = getStream(obj)
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
clear all;
warning('off','all');

% parameters & settings
k = 250; % spring stiffness
dt = 0.01; % time step
frame_num = 5000; % number of total frames to run
mass_num = 10; % number of masses in the mass-spring system

% construct a 2-dimensional mass-spring system
V = [linspace(1,1+mass_num,mass_num)',zeros(mass_num,1)];
E = [linspace(1,mass_num-1,mass_num-1)',linspace(2,mass_num,mass_num-1)']; % edges
V = V-(max(V)+min(V))/2; V = V/max(V(:)); % normalize the mesh s.t. it has the right size

vec = @(X) X(:); % vectorization: [x1x2 y1y2] order

% boundary condition (feel free to play with it!)
ci = [1 mass_num]; % fixed vertex indices
ci = [ci ci+mass_num]'; % match the [x1x2 y1y2] order --> Q: why is this important?

M = eye(2*size(V,1)); % mass matrix
g = 1*vec(repmat([0 -9.8],size(V,1),1)); % gravity

% set up the viewer
clf;
hold on;
ts = line('XData',V(:,1),'YData',V(:,2),'LineWidth',2); % draw springs as lines
tp = scatter(V(:,1),V(:,2),'.b','SizeData',400); % draw masses as dots
hold off;
axis equal;
axis([-1.5 1.5 -2 0.2]);
axis manual;
drawnow;

% initialization
P = vec(V); % mass positions for the next time step
Pt = vec(V); % mass positions for the current time step
Ptt = vec(V); % mass positions for the previous time step

for iter = 1:frame_num

  % save the states of previous time steps
  Ptt = Pt; Pt = P;

  % implicit euler time stepper
  max_iter = 50;
  for i = 1 : max_iter

    % gradient and hessian for the elastic potential energy
    [G,K] = mass_spring_gradient_hessian(V,E,k,reshape(P,size(V)));

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% solution %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % TASK: Use newton's method to solve the following nonlinear optimization problem
    %   min total_energy(P)
    %       s.t. P(ci,:) = 0
    %
    % White list:
    % - "min_quad_with_fixed" to solve the constrained optimization problem
    %
    % Hint:
    % - total energy = incremental kinetic energy + gravitational
    % potential energy + elastic potential energy

    dP = zeros(size(P));
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if norm(dP) < 1e-6
        break;
    end

    alpha = 1; % fixed step size for newton
    P = P + alpha * dP; % update P

  end

  Vnew = reshape(P,size(V));
  ts.XData = Vnew(:,1); ts.YData = Vnew(:,2); % update springs
  tp.XData = Vnew(:,1); tp.YData = Vnew(:,2); % update masses
  title(sprintf('%d : %d',iter, i),'Fontsize',20);
  drawnow;
%   figgif('./mass-spring.gif'); % uncomment this line to save a .gif

end

author: odedstein

repo: sgi-introduction-course


text:
function nearest = find_closest1 ( m, nr, r, ns, s )

%*****************************************************************************80
%
%% FIND_CLOSEST1 finds the nearest R point to each S point.
%
%  Discussion:
%
%    We are given R, a set of NR points in M dimensions.
%
%    We are given S, a set of NS points in M dimensions.
%
%    For each S(I) in S, we seek the index J of the point R(J)
%    which is nearest to S(I) over all points in R.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer NR, the number of cell generators.
%
%    Input, real R(M,NR), the cell generators.
%
%    Input, integer NS, the number of sample points.
%
%    Input, real S(M,NS), the points to be checked.
%
%    Output, integer NEAREST(NS), the index of the nearest cell generators.
%
  nearest = zeros ( ns, 1 );

  for js = 1 : ns

    distance = Inf;
    nearest(js) = -1;

    for jr = 1 : nr

      dist_sq = 0.0;
      for i = 1 : m
        dist_sq = dist_sq + ( r(i,jr) - s(i,js) )^2;
      end

      if ( dist_sq < distance )
        distance = dist_sq;
        nearest(js) = jr;
      end

    end

  end

  return
end


text:
function [d]= calc_corrt_test(d,ntrain,ntest,alpha,loss_type,as_vector)

%  CALC_corrt_test calculate corrt_test test across data objects
%
%  [X]=CALC_corrt_test(D,L,K)
%  A loss function L can also be supplied if a loss function has not
%  already been applied.
%  The (optional) parameter K if set to 0 (default) specifies if
%  the results should be stored in a data object or if set to 1, a vector.

%  The data objects D should be stored as a cell array of cell arrays
%  such that length(D) is the number of methods to compare and
%  length(D{i}) is the number of trials, which is equal for all i.
%  If data is not grouped (it is stored as a flat cell array, an attempt
%  is made to group the data automatically (e.g splitting into cv folds).

if nargin<6 as_vector=0; end;
if nargin<5 calc_loss=-1; loss_type='class_loss'; else calc_loss=1; end;
if isempty(loss_type) calc_loss=0; end;

if calc_loss==-1  %% check to see if we should calc loss
  dd=group2cell(d.child{1});
  if ~(isa(dd,'data') | isa(dd,'data_global') )  dd=dd{1}; end;
  if isempty(findstr('_loss',dd.name))
    disp('[assuming class_loss]'); calc_loss=1;
  end
end

if calc_loss==1
  d=train(loss(loss_type),d);
end

l=length(d.child);             % number of folds
dd=[];
for i=1:l
  dd{i}=group2cell(d.child{i}); %% calc elements in each fold
end
m=length(dd{1}); miscell=iscell(dd{1});

r=[]; xs=[];
for i=1:m
  r1=[];
  for j=1:m
    v=[]; v1=[]; v2=[];
    for k=[1:l]
      v1=[v1 dd{k}{i}.Y];
      v2=[v2 dd{k}{j}.Y];
      v=[v dd{k}{i}.Y-dd{k}{j}.Y];
    end
    %[p side]= corrt_test_test(v); if side==-1 p=1-p; end;
    %p=signrank(v1,v2); if mean(v1-v2)>0  p=1-p; end;
    p=correst(v1,v2,l,ntrain,ntest,alpha);
    r1=[r1 p];
  end
  xs=[xs;r1];

   % ------ remove "fold" word and add new loss value ---------

  r{i}=dd{1}{i}.name; y=r{i};
  t=findstr('fold=',y);
  if ~isempty(y) y1=[y(1:t-1) y(t+4:length(y))]; end;
  t=max(findstr('->',y));
  if ~isempty(t) %% remove loss  function calculation
    y=[y(1:t-1) '->'];
  end
  y=[y  ' corrt_test =['];
  for h=1:m
    y=[y num2str(r1(h)) ' '];
  end; y=[y ']'];
  r{i}=y;
end

d=data(char(r),xs,[]);

if as_vector
  d=d.X;
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function i4mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% I4MAT_PRINT_SOME prints out a portion of an I4MAT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of the matrix.
%
%    Input, integer A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, an optional title.
%
  incx = 10;

  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  for j2lo = max ( jlo, 1 ): incx : min ( jhi, n )

    j2hi = j2lo + incx - 1;
    j2hi = min ( j2hi, n );
    j2hi = min ( j2hi, jhi );

    inc = j2hi + 1 - j2lo;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col: ' );

    for j = j2lo : j2hi
      fprintf ( 1, '%7d  ', j );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row\n' );
    fprintf ( 1, '\n' );

    i2lo = max ( ilo, 1 );
    i2hi = min ( ihi, m );

    for i = i2lo : i2hi

      fprintf ( 1, '%5d  ', i );
      for j = j2lo : j2hi
        fprintf ( 1, '%7d  ', a(i,j) );
      end

      fprintf ( 1, '\n' );

    end

  end

  return
end


text:
%script cleanf1
if (~I(i,j))&(~I(i,j-1))...
    &I(i,j+1)&I(i-1,j)...
    &I(i-1,j+1)&(~I(i+1,j))...
    &I(i+1,j-1)&I(i+1,j+1)
  I(i,j)=1;
end


text:
function x = lprec1_old(c, d, h, g, extmod)
% LPREC1_OLD   Laplacian pyramid reconstruction using the old method
%
%	x = lprec1_old(c, d, h, g)
%
% Input:
%   c:      coarse signal at half size
%   d:      detail signal at full size
%   h, g:   two biorthogonal 1-D lowpass filters
%   extmod: [optional] extension mode (default is 'per')
%
% Output:
%   x:      reconstructed signal
%
% See also:	LPDEC1

if ~exist('extmod', 'var')
    extmod = 'per';
end

nd = ndims(c);

% Even size filter needs to be adjusted to obtain perfect reconstruction
adjust = mod(length(g) + 1, 2);

% Upsample and filter the coarse signal
p = c;
for dim = 1:nd
    p = upfilt(p, g, dim, extmod, adjust);
end

% Add with the detail signal
x = p + d;


text:
function bel = tree_protocol(engine, evidence, pot_type, local_kernel, msg)

bnet = bnet_from_engine(engine);
ns = bnet.node_sizes;
onodes = find(~isemptycell(evidence));

ndoms = length(engine.gdl.doms);
prod_of_msg = cell(1, ndoms);
bel = cell(1, ndoms);

% collect to root (node to parents)
for n=engine.postorder
  % absorb msgs from children
  prod_of_msg{n} = mk_initial_pot(pot_type, engine.gdl.doms{n}, ns, bnet.cnodes, onodes);
  for c=children(engine.tree, n)
    prod_of_msg{n} = multiply_by_pot(prod_of_msg{n}, msg{c,n});
  end
  % send msg to parents
  for p=parents(engine.tree, n)
    if iter==1
      temp = prod_of_msg{n};
    else
      temp = divide_by_pot(prod_of_msg{n}, old_msg{p,n});
    end
    temp = multiply_by_pot(temp, local_kernel{n});
    temp2 = marginalize_pot(temp, engine.gdl.sepset{n,p}, engine.maximize);
    %fprintf('%d sends %d\n', n, p);
    msg{n,p} = normalize_pot(temp2);
  end
end

% distribute from root (node to children)
for n=engine.preorder
  % absorb from parents
  %prod_of_msg{n} = mk_initial_pot(pot_type, doms{n}, ns, cnodes, onodes);
  for p=parents(engine.tree, n)
    prod_of_msg{n} = multiply_by_pot(prod_of_msg{n}, msg{p,n});
  end
  bel{n} = normalize_pot(multiply_pots(prod_of_msg{n}, local_kernel{n}));
  % send msg to children
  for c=children(engine.tree, n)
    temp = divide_by_pot(prod_of_msg{n}, msg{c,n});
    temp = multiply_by_pot(temp, local_kernel{n});
    temp2 = marginalize_pot(temp, engine.gdl.sepset{n,c}, engine.maximize);
    %fprintf('%d sends %d\n', n, c);
    msg{n,c} = normalize_pot(temp2);
  end
end


text:
function [Knots_n] = refine_linear_grid_3d(Knots, old_spacing, ds, volsz, varargin)
    %[Knots_n] = refine_linear_grid_3d(Knots, old_spacing, ds, volsz,
    %  {k_up, upsampling_type, old_size, new_spacing}
    k_down = 0.5;
    if nargin >= 5
        k_down = varargin{1};
    end
    upsampling_type = 'sample_std';
    if nargin >= 6
        upsampling_type = varargin{2};
    end
    new_spacing = old_spacing;
    if nargin >= 8
        new_spacing = varargin{4};
    end

    if false %k_down == 0.5
        Knots_sz = ceil(volsz ./ old_spacing) + 1;
        Knots_n = linear_disp_3d(Knots, Knots_sz, ds+1);
    else
        interp_type = 0;
        tmp = size(Knots);
        ksz2 = (fl(tmp(1:3)) .* (ds(:) + 1))';
        ksz2 = ceil(volsz ./ new_spacing) + 1;
    %     ksz2 = size(Knots_n); ksz2 = ksz2(1:3);

%         tmp = cat(4, volresize(Knots(:,:,:,1), ksz2, interp_type, upsampling_type), ...
%                      volresize(Knots(:,:,:,2), ksz2, interp_type, upsampling_type), ...
%                      volresize(Knots(:,:,:,3), ksz2, interp_type, upsampling_type));
        %Knots_n = tmp;

        Kt = zeros([ksz2, 3]);
        ksz3 = round([size(Knots,1),size(Knots,2),size(Knots,3)] ./ k_down);
        ksz3 = ksz3(1:3);
        tmp = cat(4, volresize_t2(Knots(:,:,:,1), 1./k_down, interp_type), ...
                     volresize_t2(Knots(:,:,:,2), 1./k_down, interp_type), ...
                     volresize_t2(Knots(:,:,:,3), 1./k_down, interp_type));
        Kt(1:size(tmp,1), 1:size(tmp,2), 1:size(tmp,3), :) = tmp;
        Kt = Kt(1:ksz2(1), 1:ksz2(2), 1:ksz2(3), :);
        Knots_n = Kt;
%         ksz2
%         ksz3
%         size(Knots_n)
%         pause;
    %     Knots_n = tmp(1:ksz1(1), 1:ksz1(2), 1:ks z1(3), :);
    end
    if strcmp(upsampling_type, 'variation')
        old_volsz = varargin{2};
        new_spacing = varargin{3};

        Tmin = linear_disp_3d(Knots, old_volsz, old_spacing);
        Tmin_u = cat(4, volresize(Tmin(:,:,:, 1),volsz, 0), volresize(Tmin(:,:,:, 2), volsz, 0), volresize(Tmin(:,:,:, 3), volsz, 0));

        objf = @(x) align_knots(x, Tmin_u, volsz, size(Knots_n), new_spacing);
        uopt = []; uopt.method = 'lbfgs'; uopt.MaxIter = 30; uopt.Corr=15; uopt.Display = 'off';
        K2 = minFunc(objf, Knots_n(:), uopt);
        Knots2 = reshape(K2, size(Knots_n));
        Knots_n = Knots2;
    end

%     imagesc([ Knots_n2(:,:, 1), Knots_n(:,:, 1)]); pause;
%     Knots_n = Knots_n * 0.9 + Knots_n2 * 0.1;
%     Knots_n = Knots_n2;
%     Knots_n = cat(4, volresize(Knots(:,:,:,1), Knots_sz, 0), ...
%         volresize(Knots(:,:,:,2), Knots_sz, 0), volresize(Knots(:,:,:,3), Knots_sz, 0));
%     Knots_n = mex_3d_linear_disp_double(Knots, Knots_sz, ds+1);

%     for i = 1 : 3
%         if all( fl(Knots_n(:,:,end, i)) == 0)
%             Knots_n(:,:,end,i) = Knots_n(:,:,end-1,i);
%         end
%     end
%     for i = 1 : 3
%         if all( fl(Knots_n(:,end,:, i)) == 0)
%             Knots_n(:,end,:,i) = Knots_n(:,end-1,:,i);
%         end
%     end
%     for i = 1 : 3
%         if all( fl(Knots_n(end,:,:, i)) == 0)
%             Knots_n(end,:,:,i) = Knots_n(end-1,:,:,i);
%         end
%     end

%     T = mex3d_linear_disp(Knots, volsz, old_spacing);
%     ksz = ceil(volsz ./ new_spacing) + 1;
% %     1+new_spacing(1):new_spacing(1):(ksz(1)-1)*new_spacing(1)
%     m1 = min((ksz(1)-1)*new_spacing(1), volsz(1));
%     m2 = min((ksz(2)-1)*new_spacing(2), volsz(2));
%     m3 = min((ksz(3)-1)*new_spacing(3), volsz(3));
%     Knots_n =...
%         T([1, round(1+new_spacing(1):new_spacing(1):m1), volsz(1)], ...
%           [1, round(1+new_spacing(2):new_spacing(2):m2), volsz(2)], ...
%           [1, round(1+new_spacing(3):new_spacing(3):m3), volsz(3)], :);

end

function [f, gr] = align_knots(K, T, szv, szk, grid_spacing)
    K = reshape(K, szk);
    Kx = linear_disp_3d(K, szv, grid_spacing);
    df = Kx - T;
    f = sum(df(:).^2)/2;
    [gr1, gr2, gr3] = linear_partial_conv_3d(df(:,:,:, 1), df(:,:,:, 2), df(:,:,:, 3), ...
                     size(K), grid_spacing);
    gr = [gr1(:); gr2(:); gr3(:)];
end

repo: pTVreg


text:
%% Learns the weights of a perceptron and displays the results.
function [w] = learn_perceptron(neg_examples_nobias,pos_examples_nobias,w_init,w_gen_feas)
%%
% Learns the weights of a perceptron for a 2-dimensional dataset and plots
% the perceptron at each iteration where an iteration is defined as one
% full pass through the data. If a generously feasible weight vector
% is provided then the visualization will also show the distance
% of the learned weight vectors to the generously feasible weight vector.
% Required Inputs:
%   neg_examples_nobias - The num_neg_examples x 2 matrix for the examples with target 0.
%       num_neg_examples is the number of examples for the negative class.
%   pos_examples_nobias - The num_pos_examples x 2 matrix for the examples with target 1.
%       num_pos_examples is the number of examples for the positive class.
%   w_init - A 3-dimensional initial weight vector. The last element is the bias.
%   w_gen_feas - A generously feasible weight vector.
% Returns:
%   w - The learned weight vector.
%%

%Bookkeeping
num_neg_examples = size(neg_examples_nobias,1);
num_pos_examples = size(pos_examples_nobias,1);
num_err_history = [];
w_dist_history = [];

%Here we add a column of ones to the examples in order to allow us to learn
%bias parameters.
neg_examples = [neg_examples_nobias,ones(num_neg_examples,1)];
pos_examples = [pos_examples_nobias,ones(num_pos_examples,1)];

%If weight vectors have not been provided, initialize them appropriately.
if (~exist('w_init','var') || isempty(w_init))
    w = randn(3,1);
else
    w = w_init;
end

if (~exist('w_gen_feas','var'))
    w_gen_feas = [];
end

%Find the data points that the perceptron has incorrectly classified
%and record the number of errors it makes.
iter = 0;
[mistakes0, mistakes1] = eval_perceptron(neg_examples,pos_examples,w);
num_errs = size(mistakes0,1) + size(mistakes1,1);
num_err_history(end+1) = num_errs;
fprintf('Number of errors in iteration %d:\t%d\n',iter,num_errs);
fprintf(['weights:\t', mat2str(w), '\n']);
plot_perceptron(neg_examples, pos_examples, mistakes0, mistakes1, num_err_history, w, w_dist_history);
key = input('<Press enter to continue, q to quit.>', 's');
if (key == 'q')
    return;
end

%If a generously feasible weight vector exists, record the distance
%to it from the initial weight vector.
if (length(w_gen_feas) ~= 0)
    w_dist_history(end+1) = norm(w - w_gen_feas);
end

%Iterate until the perceptron has correctly classified all points.
while (num_errs > 0)
    iter = iter + 1;

    %Update the weights of the perceptron.
    w = update_weights(neg_examples,pos_examples,w);

    %If a generously feasible weight vector exists, record the distance
    %to it from the current weight vector.
    if (length(w_gen_feas) ~= 0)
        w_dist_history(end+1) = norm(w - w_gen_feas);
    end

    %Find the data points that the perceptron has incorrectly classified.
    %and record the number of errors it makes.
    [mistakes0, mistakes1] = eval_perceptron(neg_examples,pos_examples,w);
    num_errs = size(mistakes0,1) + size(mistakes1,1);
    num_err_history(end+1) = num_errs;

    fprintf('Number of errors in iteration %d:\t%d\n',iter,num_errs);
    fprintf(['weights:\t', mat2str(w), '\n']);
    plot_perceptron(neg_examples, pos_examples, mistakes0, mistakes1, num_err_history, w, w_dist_history);
    key = input('<Press enter to continue, q to quit.>', 's');
    if (key == 'q')
        break;
    end
end

%WRITE THE CODE TO COMPLETE THIS FUNCTION
function [w] = update_weights(neg_examples, pos_examples, w_current)
%%
% Updates the weights of the perceptron for incorrectly classified points
% using the perceptron update algorithm. This function makes one sweep
% over the dataset.
% Inputs:
%   neg_examples - The num_neg_examples x 3 matrix for the examples with target 0.
%       num_neg_examples is the number of examples for the negative class.
%   pos_examples- The num_pos_examples x 3 matrix for the examples with target 1.
%       num_pos_examples is the number of examples for the positive class.
%   w_current - A 3-dimensional weight vector, the last element is the bias.
% Returns:
%   w - The weight vector after one pass through the dataset using the perceptron
%       learning rule.
%%
w = w_current;
num_neg_examples = size(neg_examples,1);
num_pos_examples = size(pos_examples,1);
for i=1:num_neg_examples
    this_case = neg_examples(i,:);
    x = this_case'; %Hint
    activation = this_case*w;
    if (activation >= 0)
        %YOUR CODE HERE
        w = w - x;
    end
end
for i=1:num_pos_examples
    this_case = pos_examples(i,:);
    x = this_case';
    activation = this_case*w;
    if (activation < 0)
        %YOUR CODE HERE
        w = w + x;
    end
end

function [mistakes0, mistakes1] =  eval_perceptron(neg_examples, pos_examples, w)
%%
% Evaluates the perceptron using a given weight vector. Here, evaluation
% refers to finding the data points that the perceptron incorrectly classifies.
% Inputs:
%   neg_examples - The num_neg_examples x 3 matrix for the examples with target 0.
%       num_neg_examples is the number of examples for the negative class.
%   pos_examples- The num_pos_examples x 3 matrix for the examples with target 1.
%       num_pos_examples is the number of examples for the positive class.
%   w - A 3-dimensional weight vector, the last element is the bias.
% Returns:
%   mistakes0 - A vector containing the indices of the negative examples that have been
%       incorrectly classified as positive.
%   mistakes0 - A vector containing the indices of the positive examples that have been
%       incorrectly classified as negative.
%%
num_neg_examples = size(neg_examples,1);
num_pos_examples = size(pos_examples,1);
mistakes0 = [];
mistakes1 = [];
for i=1:num_neg_examples
    x = neg_examples(i,:)';
    activation = x'*w;
    if (activation >= 0)
        mistakes0 = [mistakes0;i];
    end
end
for i=1:num_pos_examples
    x = pos_examples(i,:)';
    activation = x'*w;
    if (activation < 0)
        mistakes1 = [mistakes1;i];
    end
end

repo: neural-nets


text:
function [ x, seed ] = simplex_unit_sample ( dim_num, n, seed )

%*****************************************************************************80
%
%% SIMPLEX_UNIT_SAMPLE samples the unit simplex.
%
%  Discussion:
%
%    The interior of the unit DIM_NUM-dimensional simplex is the set of
%    points X(1:DIM_NUM) such that each X(I) is nonnegative, and
%    sum(X(1:DIM_NUM)) <= 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 August 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Reuven Rubinstein,
%    Monte Carlo Optimization, Simulation, and Sensitivity
%    of Queueing Networks,
%    Wiley, 1986.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the space.
%
%    Input, integer N, the number of points.
%
%    Input/output, integer SEED, a seed for the random number generator.
%
%    Output, real X(DIM_NUM,N), the points.
%

%
%  The construction begins by sampling DIM_NUM+1 points from the
%  exponential distribution with parameter 1.
%
  for j = 1 : n

    [ e, seed ] = r8vec_uniform_01 ( dim_num+1, seed );

    e(1:dim_num+1) = -log ( e(1:dim_num+1) );

    x(1:dim_num,j) = e(1:dim_num)' / sum ( e(1:dim_num+1) );

  end

  return
end


text:
% RES = upConv(IM, FILT, EDGES, STEP, START, STOP, RES)
%
% Upsample matrix IM, followed by convolution with matrix FILT.  These
% arguments should be 1D or 2D matrices, and IM must be larger (in
% both dimensions) than FILT.  The origin of filt
% is assumed to be floor(size(filt)/2)+1.
%
% EDGES is a string determining boundary handling:
%    'circular' - Circular convolution
%    'reflect1' - Reflect about the edge pixels
%    'reflect2' - Reflect, doubling the edge pixels
%    'repeat'   - Repeat the edge pixels
%    'zero'     - Assume values of zero outside image boundary
%    'extend'   - Reflect and invert
%    'dont-compute' - Zero output when filter overhangs OUTPUT boundaries
%
% Upsampling factors are determined by STEP (optional, default=[1 1]),
% a 2-vector [y,x].
%
% The window over which the convolution occurs is specfied by START
% (optional, default=[1,1], and STOP (optional, default =
% step .* (size(IM) + floor((start-1)./step))).
%
% RES is an optional result matrix.  The convolution result will be
% destructively added into this matrix.  If this argument is passed, the
% result matrix will not be returned. DO NOT USE THIS ARGUMENT IF
% YOU DO NOT UNDERSTAND WHAT THIS MEANS!!
%
% NOTE: this operation corresponds to multiplication of a signal
% vector by a matrix whose columns contain copies of the time-reversed
% (or space-reversed) FILT shifted by multiples of STEP.  See corrDn.m
% for the operation corresponding to the transpose of this matrix.

% Eero Simoncelli, 6/96.  revised 2/97.

function result = upConv(im,filt,edges,step,start,stop,res)

%% THIS CODE IS NOT ACTUALLY USED! (MEX FILE IS CALLED INSTEAD)

fprintf(1,'WARNING: You should compile the MEX version of "upConv.c",\n         found in the MEX subdirectory of matlabPyrTools, and put it in your matlab path.  It is MUCH faster, and provides more boundary-handling options.\n');

%------------------------------------------------------------
%% OPTIONAL ARGS:

if (exist('edges') == 1)
  if (strcmp(edges,'reflect1') ~= 1)
    warning('Using REFLECT1 edge-handling (use MEX code for other options).');
  end
end

if (exist('step') ~= 1)
  step = [1,1];
end

if (exist('start') ~= 1)
  start = [1,1];
end

% A multiple of step
if (exist('stop') ~= 1)
  stop = step .* (floor((start-ones(size(start)))./step)+size(im))
end

if ( ceil((stop(1)+1-start(1)) / step(1)) ~= size(im,1) )
  error('Bad Y result dimension');
end
if ( ceil((stop(2)+1-start(2)) / step(2)) ~= size(im,2) )
  error('Bad X result dimension');
end

if (exist('res') ~= 1)
  res = zeros(stop-start+1);
end

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

tmp = zeros(size(res));
tmp(start(1):step(1):stop(1),start(2):step(2):stop(2)) = im;

result = rconv2(tmp,filt) + res;

repo: SelfExSR


text:
function z = GetVolt()
%
%
w = 0 + 4*randn(1,1);
z = 14.4 + w;

author: philbooks

repo: Kalman-Filter-for-Beginners


text:
function [lat, lon, alt, vel, att, pos_std, vel_std, att_std, ins_status, pos_type, time] = inspvaxa2pvax(file_name)
%     data = importdata(file_name);
    fid = fopen(file_name, 'rt');

    lat = [];
    lon = [];
    alt = [];
    vel = [];
    att = [];
    pos_std = [];
    vel_std = [];
    att_std = [];
    ins_status = [];
    pos_type = [];
    time = [];

%     for i = 1:length(data)
%         str = char(data(i));
    while ~feof(fid)
        str = fgetl(fid);
        if length(str)>10 && (strcmp(str(1:9),'#INSPVAXA'))
            sstr = string(str);
            sstr = sstr.split(';');

            sstr1 = sstr(1);
            sstr1 = sstr1.split(',');
            gps_week_time = mod(str2double(sstr1(7)), (3600*24)) - 18;
            time = [time; gps_week_time];

            sstr2 = sstr(2);
            sstr2 = sstr2.split(',');

            lat = [lat; str2double(sstr2(3))];
            lon = [lon; str2double(sstr2(4))];
            alt = [alt; str2double(sstr2(5))];
            vel = [vel; str2double(sstr2([8,7,9]))'];
            att = [att; str2double(sstr2([11,10,12]))'];
            pos_std = [pos_std; str2double(sstr2(13:15))'];
            vel_std = [vel_std; str2double(sstr2([17,16,18]))'];
            att_std = [att_std; str2double(sstr2([20,19,21]))'];

            ins_status = [ins_status; sstr2(1)];
            pos_type = [pos_type; sstr2(2)];
        end
    end

    fclose(fid);
end

author: yandld


text:
classdef TestCalcOpticalFlowDF
    %TestCalcOpticalFlowDF

    methods (Static)
        function test_1
            im1 = 255*uint8([...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 1 1 1 0 0 0 0;...
                0 0 0 1 0 1 0 0 0 0;...
                0 0 0 1 1 1 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
                0 0 0 0 0 0 0 0 0 0;...
            ]);
            im2 = circshift(im1, [0 1]);
            flow = cv.calcOpticalFlowDF(im1, im2);
            validateattributes(flow, {'single'}, ...
                {'3d', 'size',[size(im1,1) size(im1,2) 2]});
        end

        function test_2
            prevImg = cv.imread(fullfile(mexopencv.root(),'test','RubberWhale1.png'), ...
                'Grayscale',true, 'ReduceScale',2);
            nextImg = cv.imread(fullfile(mexopencv.root(),'test','RubberWhale2.png'), ...
                'Grayscale',true, 'ReduceScale',2);
            flow = cv.calcOpticalFlowDF(prevImg, nextImg);
            validateattributes(flow, {'single'}, ...
                {'3d', 'size',[size(prevImg,1) size(prevImg,2) 2]});
        end

        function test_error_argnum
            try
                cv.calcOpticalFlowDF();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%
% Optimize the order of strokes in a MotorProgram,
% before the relations are fixed
%
% Input
%  Q: reference to MotorProgram object
%
function optimize_order_MP(Q,lib)

    ps = defaultps_bottomup;

    % get all the permutations to try
    if Q.ns <= ps.max_ns_all_perm
        % try all possible combinations
        P = perms(1:Q.ns);
    else
        % try a subset of the permutations
        np = factorial(ps.max_ns_all_perm); %720
        P = zeros(np,Q.ns);
        for i=1:np
            P(i,:) = randperm(Q.ns);
        end
        P = unique(P,'rows');
    end

    % score all the permutations
    n = size(P,1);
    scores = zeros(n,1);
    for i=1:n
       QQ = Q.copy();
       perm = P(i,:);
       QQ.S = QQ.S(perm);
       scores(i) =  scoreMP_NoRel(QQ,lib,'type',true,'token',true,'stat',false,'image',false);
    end

    % pick the best order, and return Q without instantiating relations
    [~,windx] = randmax(scores);
    perm = P(windx,:);
    Q.S = Q.S(perm);
end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
conf = config;

load([conf.root_path, 'splits.mat'], 'train_images', 'val_images');

% check if all images are readible
%check_images(conf.image_path, train_images);

training_images = load_cropped_images(conf.image_path, train_images, ...
    conf.preprocessing.crop_padding, conf.preprocessing.square_size);

savefast(fullfile(conf.root_path,...
    'training_images_crop15_square256.mat'), 'training_images');

validation_images = load_cropped_images(conf.image_path, val_images, ...
    conf.preprocessing.crop_padding, conf.preprocessing.square_size);

savefast(fullfile(conf.root_path,...
    'validation_images_crop15_square256.mat'), 'validation_images');

author: rksltnl


text:
function [matlabclass,bytes_per_element] = dtypechar2matlabclass(dtypechar)

switch dtypechar,
  case {'c','s','p'},
    matlabclass = 'char';
    bytes_per_element = 1;
  case 'b',
    matlabclass = 'int8';
    bytes_per_element = 1;
  case 'B',
    matlabclass = 'uint8';
    bytes_per_element = 1;
  case 'h',
    matlabclass = 'int16';
    bytes_per_element = 2;
  case 'H',
    matlabclass = 'uint16';
    bytes_per_element = 2;
  case {'i','l'},
    matlabclass = 'int32';
    bytes_per_element = 4;
  case {'I','L'},
    matlabclass = 'uint32';
    bytes_per_element = 4;
  case 'q',
    matlabclass = 'int64';
    bytes_per_element = 8;
  case 'Q',
    matlabclass = 'uint64';
    bytes_per_element = 8;
  case 'f',
    matlabclass = 'float';
    bytes_per_element = 4;
  case 'd',
    matlabclass = 'double';
    bytes_per_element = 8;
  otherwise
    error('Unknown data type %s',dtypechar);
end

author: kristinbranson

repo: JAABA


text:
function n=im3Dnorm(img,normind,varargin)
% IMAGE3DNORM computes the desired image norm
%   IMAGE3DNORM(IMG,NORMIND) computes the norm if image IMG using the norm
%   defined in NORMING
%
%   IMG         A 3D image
%   NORMIND     'LN': N norm
%               'TV': TV norm
%
%
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri
%--------------------------------------------------------------------------
if ~ischar(normind)
      error('CBCT:image3Dnorm:InvalidInput','Norm option has to be a string');
end
if length(normind)~=2
      error('CBCT:image3Dnorm:InvalidInput','Unknown norm option');
end

if strcmp(normind(1),'L')
    Nnorm=str2double(normind(2:end));
    n=norm(img(:),Nnorm);
    return;
end

if strcmp(normind,'TV')
    typeGrad=varargin{1};
    if strcmpi(typeGrad,'central')
        [gx,gy,gz]=gradient(img);
    end
    if strcmpi(typeGrad,'forward')
        gx=diff(img,1,1);
        gy=diff(img,1,2);
        gz=diff(img,1,3);
        gx=cat(1,gx,zeros(size(gx(end,:,:))));
        gy=cat(2,gy,zeros(size(gy(:,end,:))));
        gz=cat(3,gz,zeros(size(gz(:,:,end))));

    end
    if strcmpi(typeGrad,'backward')
        gx=diff(img,1,1);
        gy=diff(img,1,2);
        gz=diff(img,1,3);
        gx=cat(1,zeros(size(gx(1,:,:))),gx);
        gy=cat(2,zeros(size(gy(:,1,:))),gy);
        gz=cat(3,zeros(size(gz(:,:,1))),gz);
    end
    if ~exist('gx','var')
        error('CBCT:image3Dnorm:InvalidInput','Unknown gradient option for the TV norm');
    end
    g=sqrt(gx.^2+gy.^2+gz.^2); clear gx gy gz;
    n=sum(g(:));
    return;
end

error('CBCT:image3Dnorm:InvalidInput','Unknown norm option');
end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function gZ = vrate(Z,s)
% PURPOSE: Compute the year-on-year rate of growth of a vector time series
% ------------------------------------------------------------------------
% SYNTAX:  gZ = vrate(Z,s);
% ------------------------------------------------------------
% OUTPUT: gZ: nxk --> yoy rate of growth. First s obs are NaN
% ------------------------------------------------------------
% INPUT:  Z: nxk --> vector time series to be filtered
%         s: 1x1 --> number of periods per year
% ------------------------------------------------------------
% LIBRARY: rate
% -----------------------------------------------------------------------

% written by:
%  Enrique M. Quilis
%  Macroeconomic Research Department
%  Ministry of Economy and Competitiveness
%  <enrique.quilis@mineco.es>

% Version 1.1 [August 2006]

[n,k] = size(Z);

% Preallocation
gZ = NaN * ones(n,k);

% Basic loop
for j=1:k
    gz = rate(Z(:,j),s);
    gZ(:,j) = gz;
end

repo: mcbench-benchmarks


text:
function tol = cputime_tol
    % cputime_res applies to either t1 or t2 individually;
    % tol applies to t2-t1;
    tol = 2*cputime_res + sqrt(eps);
end

%!test
%! tol = cputime_tol;


text: %% S2AxisField

author: mtex-toolbox

repo: mtex

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


text:
 function [kp, iop, rfp] = example_spsp_param1()
%function [kp, iop, rfp] = example_spsp_param1()
% setup parameters for spectral-spatial pulse design.
%
%Chun-yu Yip, 4/7/09

kp = kparametersSPSP();                 % k-space trajectory parameters
iop = ioparametersSPSP();               % input/ouput parameters
rfp = rfparametersSPSP();               % RF pulse design parameters

function [iop] = ioparametersSPSP()
%function [iop] = ioparametersSPSP()
%iop is a structure that contains parameters for input and
%output of the pulse design code.

iop.writetofile_sim = logical(1);       %Write designed pulse and gz waveforms
                                        %to files for Bloch simulation.

iop.writetofile_scanner = logical(1);   %Write designed pulse and gz waveforms
                                        %to files for MRI scanner.

%Bloch simulation parameters
iop.bfovz = 3;                          %cm; range of z (space) in Bloch
                                        %simulation = [-bfovz/2, bfovz/2).
iop.bdimz = 200;                        %Number of z samples in Bloch simulation
iop.bfovf = 500;                        %Hz; range of f (frequency) in Bloch
                                        %simulation = [-bfovf/2, bfovf/2).
iop.bdimf = 200;                        %Number of f samples in Bloch simulation

%File names
iop.infname = 'spspexample';            %File name for RF and z gradient waveform
                                        %files. Extension should not be included.
iop.outfname_from_sim = 'mspsp';        %File name for Bloch simulation results.
                                        %Extension should not be included.

iop.show_blochsimresults = logical(1);  %Display Bloch simulation results

iop.waveformfilespath = '.';            %Where pulse gradient waveforms and
                                        %simulation results reside.

 function [rfp] = rfparametersSPSP()
%function [rfp] = rfparametersSPSP()
%rfp is a structure that contains SPSP RF waveform computation parameters.
% Note: default parameter values as on our paper. They will generate
% the alpha=-2E-4 pulse, as shown in Fig. 3b.

    %Slice selection parameters
    rfp.slthickz = .5 ;              %cm; slice thickness. Full width
                                     %half max (FWHM) when shape is
                                     %"gaussian".
    rfp.flipangle = 30;              %degrees; flip angle.
    rfp.profileshape = 'rect';       %'rect' (recommended) or 'gaussian'
    rfp.smoothprofile = logical(0);  %Smoothing of the slice profile applied
                                     %to the desired pattern (not recommended).
    rfp.kernelstdv = 0.3;            %cm; Gaussian smoothing kernel std dev.
                                     %This parameter is useless when
                                     %rfp.smoothprofile = logical(0).
    rfp.hamming = logical(1);        %Hamming window in k-space applied
                                     %directly to the designed pulse, as
                                     %described in paper. It is intended to
                                     %smooth the excited slice profile.
                                     %Recommended, although theoretically
                                     %suboptimal.

    %Desired pattern parameters
    rfp.dfovz = 20;                  %cm; range of desired pattern definition
                                     %along z (space). d defined over [-dfovz/2,
                                     %dfovz/2).
    rfp.ddimz = 1200;                %Number of samples along z.
    rfp.dfovf = 500;                 %Hz; range of desired pattern definition
                                     %along f (frequency). f defined over
                                     %[-dfovf/2, dfovf/2).
    rfp.ddimf = 300;                 %Number of samples along f.

    %Phase pattern parameters
    rfp.zshift = -rfp.slthickz/2;    %cm; z-direction shift applied to phase
                                     %pattern. Refer to Eq. 5 in paper.
    rfp.alpha = -2E-4;               %g/cm/Hz; key proportionality constant
                                     %between field offset and field gradient
                                     %Refer to Eq. 3 in paper.
    rfp.TE = 30E-3;                  %s; echo time of GRE sequence for fMRI.
                                     %Note that it is defined as from SPSP
                                     %pulse *center* to the time when the DC
                                     %sample is acquired.

    %Conjugate gradient parameter
    rfp.Niter = 100;                 %Number of CG iterations in design process
    rfp.beta = 0;                    %Tikhonov regularization parameter; controls
                                     %pulse energy.

%function [kp] = kparametersSPSP() kp is a structure that contains
%excitation k-space trajectory parameters for spectral-spatial pulse
%design. The trajectory comes from an oscillatory Gz waveform comprised of
%trapazoids with opposite polarities, which is most commonly used for
%conventional SPSP pulses.
%
% Note: default parameter values as on our paper. They will generate
% the alpha=-2E-4 pulse, as shown in Fig. 3b.

function [kp] = kparametersSPSP()
kp.pointtime = 4E-6;            %s; RF pulse sampling period,
                                %assuming it is same for z gradient.
                                %For example, for Siemens Trio, it is
                                %2.5E-6 s. For GE Signa (3T) in Michigan,
                                %it was 4E-6 s.

kp.gmax = 4;                    %g/cm; maximum gradient amplitude
kp.dgdtmax = 15000;             %g/cm/s; maximum gradient slew rate

kp.T = 0.003;                   %s; z gradient oscillation period
kp.Ntraps = 7;                  %Number of trapezoids in z gradient waveform
                                %Note:
                                %These two parameters determine the shape
                                %of the Gz waveform. Suitable values of
                                %these two values depend on 1. slice
                                %profile shape and thickness, 2. "alpha", 3.
                                %TE(and TD), and perhaps more. As I have
                                %not come up with formuli for them,
                                %I always find good values for them by trial
                                %and error (ie, put in some reasonable
                                %values, do the pulse design, and perform Bloch
                                %simulation to see if those values give me
                                %a good excitation pattern). For a certain
                                %excitation accuracy, those values should
                                %be chosen so that the pulse is as short as
                                %possible. For pulses that we show in Fig. 3
                                %of our paper, we used:
                                %   alpha=-2:   T = 0.003;    Ntraps = 7;
                                %   alpha=-2.5: T = 0.003;    Ntraps = 8;
                                %   alpha=-3:   T = 0.003;    Ntraps = 9;
                                %Other good values for different alphas:
                                %   alpha=-1:   T = 0.0025;   Ntraps = 6;
                                %   alpha=-1.25:T = 0.0025;   Ntraps = 6;
                                %   alpha=-1.5: T = 0.0025;   Ntraps = 6;
                                %   alpha=-1.75:T = 0.0025;   Ntraps = 7;
                                %   alpha=-2.25:T = 0.003;    Ntraps = 7;
                                %   alpha=-2.75:T = 0.003;    Ntraps = 8;
                                %   alpha=-3.25:T = 0.003;    Ntraps = 10;
                                %   alpha=-3.5: T = 0.003;    Ntraps = 10;

kp.pw = 9999;                   %s; pulse duration. Will be filled once the
                                %pulse is designed (9999: not known yet)
kp.npnts = 9999;                %Number of pulse samples. Will be filled once
                                %the pulse is designed (9999: not known yet)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function F = diff(F, k, dim)
%DIFF   Componentwise partial derivative of a CHEBFUN3V object.
%   DIFF(F) is the derivative of components of F along the first variable.
%
%   DIFF(F, K) is the Kth derivative of each component of F along the first
%   variabel.
%
%   DIFF(F, K, DIM) is the Kth derivative of F along the dimension DIM.
%   DIM = 1 (default) is the derivative in the 1st input variable.
%   DIM = 2 is the derivative in the 2nd variable.
%   DIM = 3 is the derivative in the 3rd varialbe.
%
%   DIFF(F, [K1 K2], [DIM1 DIM2]) means K1-th derivative of F in dimension
%   DIM1 and K2-th derivative in dimension DIM2. DIM1 and DIM2 can be 1, 2
%   or 3 in any order.
%
%   DIFF(F, [K1 K2 K3]) is the K1-th partial derivative of F in the first
%   variable, K2-th partial derivative of F in the second variable and
%   K3-th partial derivative of F in the third variable.
%   For example, DIFF(F, [1 2 3]) is d^6F/(dx d^2y d^3z).

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

% Empty check:
if ( isempty(F) )
    return
end

% Defaults:
if ( ( nargin == 1 ) || isempty(k) )
    k = 1;
end
if ( nargin < 3 )
    dim = 1;
end

% Diff each component.
for j = 1:F.nComponents
    F.components{j} = diff(F.components{j}, k, dim);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function fpalm_render(handles, h)
    outpath = get(handles.out_dir_edit, 'String'); %load output path from the edit box
    if outpath(length(outpath))~='\'
        outpath = [outpath,'\'];
    end

    mat_file = getappdata(0,'temp_runtime');

    warning off MATLAB:load:variableNotFound    %No need to warn if one isnt found, this is taken care of below
                                                %In addition, depending on
                                                %what calls this function, some of these variables may not yet exist anyway
    load(mat_file,'xf_all','yf_all','a0_all','total_molecules','file_info',...
                  'w','preview','plot1','x_size','y_size','iprod_thresh',...
                  'threshold','upper_threshold','n_bright_pixel_threshold', ...
                  'max_pixels_above_threshold','base_name','rbox');
    warning on MATLAB:load:variableNotFound     %Return warning to its original state

%---------------
%Added for backward compatability (changes in naming mechanism)
    if exist('file_info','var')
        if ~isfield(file_info,'part_name')
            load(mat_file,'base_name');
            file_info.part_name = base_name;
        end

        if ~isfield(file_info,'start')
            load(mat_file,'n_start');
            file_info.start = n_start;
        end

        if ~isfield(file_info,'stop')
            load(mat_file,'n_end');
            file_info.stop = n_end;
        end
    else
            load(mat_file,'base_name','n_end','n_start');
            file_info.part_name = base_name;
            file_info.start = n_start;
            file_info.stop = n_end;
    end
%----------------

    if isappdata(0,'from_mat_file')   %If not being called by einzelreader.m, create new waitbar instead of updating
        preview = 0;
        w = waitbarxmod(0,'Executing "fpalm render einzelreader.m" ...','CreateCancelBtn','delete(gcf)');
        set(w,'Name','Progress Bar');
        uicontrol('Style','pushbutton','Parent',w,'String','Pause','Position',[210,10,60,23], ...
                  'UserData',1,'Callback',@pause_gui);
        pause(.1); %Pause to ensure window completes drawing
        drawnow;    %Draw the extra button immediately
        keepontop('Progress Bar');
    else
        waitbarxmod(0,w,'Executing "fpalm render einzelreader.m" ...');   %rename the waitbar
    end

    if ~exist('a0_all','var') && preview  %If einzelreader.m hasn't found anything, return with no calculations
        return
    elseif ~exist('a0_all','var') && ~preview
        delete(w);
        error('No molecules found.');
    end

    sb     = get(handles.sb,'Value');
    sbloc  = str2num(get(handles.sb_loc,'String')); %str2num must be used here, nonscalar
    x0_bar = sbloc(1); %coordinates of scale bar
    y0_bar = sbloc(2);

    q      = str2double(get(handles.cam_pix_edit,'String'));  % pixel size in microns
    exf    = str2double(get(handles.exf_edit,'String')); %expansion factor
    bkgn   = str2double(get(handles.bkgn_noise,'String')); % background noise in photons
    ppp    = str2double(get(handles.pix_to_pho,'String'));

    NA     = str2double(get(handles.NA,'String'));
    wvlnth = str2double(get(handles.wvlnth,'String'))/1000; %Conversion to um
    psf_scale = str2double(get(handles.psf_scale_edit,'String'));

    r0 = psf_scale*0.55*wvlnth/NA/1.17; % 1/e2 radius of PSF, use 1/e2 = FWHM/1.17 from Pawley
                                  % scaled (20% for "real objective" due to
                                  % measured PSF from Hess and Webb 2002)
    psf_std=r0/2; %standard deviation of psf

    ishif=str2double(get(handles.ishift,'String'));
    jshif=str2double(get(handles.jshift,'String'));
    xw_render=str2double(get(handles.xw_render_edit,'String'));
    yw_render=str2double(get(handles.yw_render_edit,'String'));

    if xw_render ~= 0
        xw=xw_render;
    else
        xw=x_size;
    end
    if yw_render ~= 0
        yw=yw_render;
    else
        yw=y_size;
    end

    r0pix=r0/q; %1/e^2 radius of psf in pixels
    psf_w02=r0pix*r0pix;

    [xpix,ypix] = meshgrid(-rbox:rbox,-rbox:rbox);
    yfit=exp(-2*((xpix).*(xpix)+(ypix).*(ypix))/psf_w02);
    npix=sum(sum(yfit));     % area of molecule in square pixels

    a0_phot=a0_all/ppp;    % peak amplitude of each molecule in photons
    N=npix*a0_phot;        % number of photons for each molecule

    lp2=((psf_std^2)+(q^2)/12)*1./N+8*pi*(psf_std^4)*(bkgn^2)/(q^2)*1./(N.*N);
    lp=sqrt(lp2);
    lppix=lp/q*exf;   % width of each molecule in pixels, based on localization-uncertainty
    lp2pix=lppix.*lppix;

    xf=exf*(xf_all-ishif-0.5);
    yf=exf*(yf_all-jshif-0.5);
    xc=uint16(xf);
    yc=uint16(yf);

    nstart=1;
    nend=total_molecules;
    render_func;

    maxintens=max(max(impts));
    impts = impts/maxintens;

    if sb %draw scale bar if selected
        um1=1/q*exf;
        um1_round=round(um1);
        bar_width=12;

        for i=1:um1_round
            for j=1:bar_width
                xi=i+x0_bar;
                yi=j+y0_bar;
                impts(yi,xi)=1;
            end
        end

        nm250=0.25/q*exf;
        nm250_round=round(nm250);
        bar_width=4;
        y0_bar=y0_bar+bar_width;

        for i=1:um1_round
            for j=1:bar_width
                ival=round(i/nm250_round);
                ival=mod(ival,2);

                xi=i+x0_bar;
                yi=j+y0_bar;
                impts(yi,xi)=double(ival);
            end
        end
    end

    if preview
        set(0,'CurrentFigure',h);
        set(h,'CurrentAxes',plot1);
        imagesc(impts,'HitTest','off') %to look at pts only
        colormap gray
        axis image;
        title(['FPALM Image: Total Molecules:',num2str(total_molecules)])
        set(plot1, 'ButtonDownFcn', {@focus_swap,guidata(gcbo)});
    else
        figure('Name','Final Image')
        imagesc(impts)
        colormap gray
        axis image;
        title(['FPALM Image: Total Molecules:',num2str(total_molecules)])
        delete(w);
    end

    if isappdata(0,'Save')
        if getappdata(0,'Save')
            if (file_info.part_name(length(file_info.part_name)) == '_') %Avoid saving with repeat seperators
                if get(handles.adv_cb,'Value')
                    impts_gray_file=[outpath,'\',file_info.part_name,'pts_',num2str(file_info.start),'-',num2str(file_info.stop),'_t',num2str(iprod_thresh),'-',num2str(threshold),'-',num2str(upper_threshold),'_npt',num2str(n_bright_pixel_threshold),'-',num2str(max_pixels_above_threshold),'_e',num2str(exf),'.tif'];
                else
                    impts_gray_file=[outpath,'\',file_info.part_name,'pts_',num2str(file_info.start),'-',num2str(file_info.stop),'_t',num2str(iprod_thresh),'_e',num2str(exf),'.tif'];
                end
            else
                if get(handles.adv_cb,'Value')
                    impts_gray_file=[outpath,'\',file_info.part_name,'_pts_',num2str(file_info.start),'-',num2str(file_info.stop),'_t',num2str(iprod_thresh),'-',num2str(threshold),'-',num2str(upper_threshold),'_npt',num2str(n_bright_pixel_threshold),'-',num2str(max_pixels_above_threshold),'_e',num2str(exf),'.tif'];
                else
                    impts_gray_file=[outpath,'\',file_info.part_name,'_pts_',num2str(file_info.start),'-',num2str(file_info.stop),'_t',num2str(iprod_thresh),'_e',num2str(exf),'.tif'];
                end
            end
            imwrite(uint8(impts*255),impts_gray_file,'Compression','none');
            rmappdata(0,'Save');

            if get(handles.oc_wf,'Value') %save widefield image if selected
                load(mat_file, 'image_sum');
                if (file_info.part_name(length(file_info.part_name)) == '_') %Avoid saving with repeated seperators
                    wf_file = [outpath,'\',file_info.part_name,'wf.tif'];
                else
                    wf_file = [outpath,'\',file_info.part_name,'_wf.tif'];
                end
                imwrite(uint8((image_sum/max(image_sum(:)))*255),wf_file,'Compression','none');
            end

            if get(handles.oc_hist,'Value') %show photon histogram of complete image
                lp_min = min(lp)*1000;
                lp_max = max(lp)*1000;
                n_bins=round(2*total_molecules^(1/3));
                bins=linspace(lp_min,lp_max,n_bins);

                figure
                hist(lp*1000,bins);
                ylabel('Frequency');
                xlabel('Localization Precision (nm)');
            end
        end
    end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB