text:
function res = ssc_W(Y,x,ta,sc,type,opC,rl)
% PURPOSE: Temporal disaggregation using the dynamic Chow-Lin method
%          proposed by Santos Silva-Cardoso (2001).
%          Without pretesting for intercept
% ------------------------------------------------------------
% SYNTAX: res = ssc_W(Y,x,ta,sc,type,opC,rl);
% ------------------------------------------------------------
% OUTPUT: res: a structure
%           res.meth    ='Santos Silva-Cardoso';
%           res.ta      = type of disaggregation
%           res.type    = method of estimation
%           res.opC     = option related to intercept
%           res.N       = nobs. of low frequency data
%           res.n       = nobs. of high-frequency data
%           res.pred    = number of extrapolations
%           res.sc      = frequency conversion between low and high freq.
%           res.p       = number of regressors (including intercept)
%           res.Y       = low frequency data
%           res.x       = high frequency indicators
%           res.y       = high frequency estimate
%           res.y_dt    = high frequency estimate: standard deviation
%           res.y_lo    = high frequency estimate: sd - sigma
%           res.y_up    = high frequency estimate: sd + sigma
%           res.u       = high frequency residuals
%           res.U       = low frequency residuals
%           res.gamma    = estimated model parameters (including y(0))
%           res.gamma_sd = estimated model parameters: standard deviation
%           res.gamma_t  = estimated model parameters: t ratios
%           res.phi      = dynamic parameter phi
%           res.beta     = estimated model parameters (excluding y(0))
%           res.beta_sd = estimated model parameters: standard deviation
%           res.beta_t  = estimated model parameters: t ratios
%           res.phi     = innovational parameter
%           res.aic     = Information criterion: AIC
%           res.bic     = Information criterion: BIC
%           res.val     = Objective function used by the estimation method
%           res.wls     = Weighted least squares as a function of phi
%           res.loglik  = Log likelihood as a function of phi
%           res.r       = grid of innovational parameters used by the estimation method
%           res.et       = elapsed time
% ------------------------------------------------------------
% INPUT: Y: Nx1 ---> vector of low frequency data
%        x: nxp ---> matrix of high frequency indicators (without intercept)
%        ta: type of disaggregation
%            ta=1 ---> sum (flow)
%            ta=2 ---> average (index)
%            ta=3 ---> last element (stock) ---> interpolation
%            ta=4 ---> first element (stock) ---> interpolation
%        sc: number of high frequency data points for each low frequency data points
%            sc= 4 ---> annual to quarterly
%            sc=12 ---> annual to monthly
%            sc= 3 ---> quarterly to monthly
%        type: estimation method:
%            type=0 ---> weighted least squares
%            type=1 ---> maximum likelihood
%        opC: 1x1 option related to intercept
%            opc = 0 : no intercept in hf model
%            opc = 1 : intercept in hf model
%        rl: innovational parameter
%            rl = []: 0x0 ---> rl=[0.05 0.99], 100 points grid
%            rl: 1x1 ---> fixed value of phi parameter
%            rl: 1x2 ---> [r_min r_max] search is performed
%                on this range, using a 200 points grid
% ------------------------------------------------------------
% LIBRARY: aggreg, criterion_SSC
% ------------------------------------------------------------
% SEE ALSO: chowlin, litterman, fernandez, td_plot, td_print
% ------------------------------------------------------------
% REFERENCE: Santos, J.M.C. and Cardoso, F.(2001) "The Chow-Lin method
% using dynamic models",Economic Modelling, vol. 18, p. 269-280.
% Di Fonzo, T. (2002) "Temporal disaggregation of economic time series:
% towards a dynamic extension", Dipartimento di Scienze Statistiche,
% Universita di Padova, Working Paper n. 2002-17.

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

% Version 2.1 [April 2012]

t0=clock;

% ------------------------------------------------------------
% Size of the problem

[N,M] = size(Y);    % Size of low-frequency input
[n,p] = size(x);    % Size of p high-frequency inputs (without intercept)

% ------------------------------------------------------------
% Initial checks

if ((opC ~= 0) & (opC ~= 1))
      error ('*** c has an invalid value ***');
end

% ------------------------------------------------------------
% Preparing the X matrix: including an intercept if opC==1

if (opC == 1)
   e=ones(n,1);
   x=[e x];       % Expanding the regressor matrix
   p=p+1;         % Number of p high-frequency inputs (plus intercept)
end

% ------------------------------------------------------------
% Generating the aggregation matrix

C = aggreg(ta,N,sc);

% -----------------------------------------------------------
% Expanding the aggregation matrix to perform
% extrapolation if needed.

if (n > sc * N)
   pred = n - sc*N;           % Number of required extrapolations
   C = [C zeros(N,pred)];
else
   pred = 0;
end

% -----------------------------------------------------------
% Temporal aggregation of the indicators

X = C*x;

% -----------------------------------------------------------
% Optimization of objective funcion (Lik. or WLS)

nrl = length(rl);

switch nrl
    case 0 %Default
        rl = [0.05 0.99];
        rex = criterion_SSC(Y,x,ta,type,opC,rl,X,C,N,n,p);
        phi = rex.phi;
        r   = rex.r;
    case 1 %Fixed innovational parameter
        phi = rl;
        type = 4;
    case 2
        rex = criterion_SSC(Y,x,ta,type,opC,rl,X,C,N,n,p);
        phi = rex.phi;
        r   = rex.r;
    otherwise
        error('*** Grid search on phi is improperly defined. Check rl ***');
end

% -----------------------------------------------------------
% Final estimation with optimal phi
% -----------------------------------------------------------

% Auxiliary matrix useful to simplify computations

I=eye(n); w=I;
LL = diag(-ones(n-1,1),-1);

% -----------------------------------------------------------
% Final estimation with optimal phi

% -----------------------------------------------------------
% Generation of difference matrix D_phi

D_phi = I + phi*LL;
D_phi(1,1)=sqrt(1-phi^2);
iD_phi = inv(D_phi);

% -----------------------------------------------------------
% Truncation remainder: q parameter

q=zeros(n,1);
q(1)=phi;

% -----------------------------------------------------------
% Expanded set of regressors: high and low frequency

z = [x q];
z_phi = iD_phi * z;
Z_phi = C * z_phi;

% -----------------------------------------------------------
% GLS estimator of gamma

w = inv(D_phi' * D_phi);
W = C * w * C';
iW = inv(W);

gamma = (Z_phi' * iW * Z_phi) \ (Z_phi' * iW * Y); % gamma GLS

U = Y - Z_phi * gamma;           % Low frequency residuals
scp = U' * iW * U;               % Weighted least squares
sigma_a = scp/(N-p-1);           % sigma_a estimator (p+1 due to lagged endogenous)
L = w * C' * iW;                 % Filtering matrix
u = L * U;                       % High frequency residuals

% -----------------------------------------------------------
% Temporally disaggregated time series

y = z_phi * gamma + u;

% -----------------------------------------------------------
% Information criteria
% p is expanded to include lagged endogenous

aic=log(sigma_a)+2*(p+1)/N;
bic=log(sigma_a)+log(N)*(p+1)/N;

% -----------------------------------------------------------
% VCV matrix of high frequency estimates

sigma_gamma = sigma_a * inv(Z_phi' * iW * Z_phi);

VCV_y = sigma_a * (eye(n)-L*C) * w + (z_phi-L*Z_phi) * sigma_gamma * (z_phi-L*Z_phi)';

d_y=sqrt((diag(VCV_y)));   % Std. dev. of high frequency estimates
y_li=y-d_y;                % Lower lim. of high frequency estimates
y_ls=y+d_y;                % Upper lim. of high frequency estimates

% -----------------------------------------------------------
% -----------------------------------------------------------
% Loading the structure

res.meth='Santos Silva-Cardoso';

% -----------------------------------------------------------
% Basic parameters

res.ta        = ta;
res.type      = type;
res.N         = N;
res.n         = n;
res.pred      = pred;
res.sc        = sc;
res.p         = p;
res.opC       = opC;
res.rl        = rl;

% -----------------------------------------------------------
% Series

res.Y         = Y;
res.x         = x;
res.y         = y;
res.y_dt      = d_y;
res.y_lo      = y_li;
res.y_up      = y_ls;

% -----------------------------------------------------------
% Residuals

res.u         = u;
res.U         = U;

% -----------------------------------------------------------
% Parameters

res.gamma     = gamma;
res.gamma_sd  = sqrt(diag(sigma_gamma));
res.gamma_t   = gamma./sqrt(diag(sigma_gamma));
res.phi       = phi;

res.beta      = gamma(1:end-1);
res.beta_sd   = res.gamma_sd(1:end-1);
res.beta_t    = res.gamma_t(1:end-1);
res.sigma_a   = sigma_a;

% -----------------------------------------------------------
% Information criteria

res.aic       = aic;
res.bic       = bic;

% -----------------------------------------------------------
% Objective function

if (nrl == 1)
   res.val      = [];
   res.wls      = [];
   res.loglik   = [];
   res.r        = [];
else
   res.val       = rex.val;
   res.wls  	 = rex.wls;
   res.loglik    = rex.loglik;
   res.r         = r;
end

% -----------------------------------------------------------
% Elapsed time

res.et  = etime(clock,t0);

author: Sable

repo: mcbench-benchmarks


text:
function next = next_boundary_node_ql ( node )

%*****************************************************************************80
%
%% NEXT_BOUNDARY_NODE_QL returns the next boundary node in a QL element.
%
%  Reference Element QL:
%
%    |
%    1  4---5---6
%    |  |       |
%    |  |       |
%    S  |       |
%    |  |       |
%    |  |       |
%    0  1---2---3
%    |
%    +--0---R---1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE, the index of the current node.  An input
%    value of 0 (or any "unusual" value") indicates that the
%    first edge node is desired.
%
%    Output, integer NEXT, the index of the next edge node.
%
  if ( node == 1 )
    next = 2;
  elseif ( node == 2 )
    next = 3;
  elseif ( node == 3 )
    next = 6;
  elseif ( node == 6 )
    next = 5;
  elseif ( node == 5 )
    next = 4;
  else
    next = 1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function rotMatrix = fick2rotMatrix(vfick,varargin)

% FICK2ROTMATRIX  transforms Fick coordinates into a 3*3 rotation matrix
%
%     ROTMATRIX = FICK2ROTMATRIX(VFICK) returns the rotation matrix ROTMATRIX
%     corresponding to the Fick coordinates VFICK [deg].
%        - VFICK is a 3-vector or a 3*N array (column i represents
%           rotation i) where N is the number of rotations. The 3-D
%           vector VFICK (or VFICK(:,i)) can be written as VFICK = [H V T]
%           where (from the subject's perspective):
%            - H is the horizontal fick angle (positive when to the left)
%            - V is the vertical fick angle (around the inter-aural axis).
%               It is positive when down.
%            - T is the torsional fick angle (around the line of signt). It
%               is positive when clockwise.
%        - ROTMATRIX is a 3*3 array or a 3*3*N tensor (in which case,
%            ROTMATRIX(:,:,ii) is a 3*3 matrix representing the rotation
%            matrix corresponding to rotation ii.
%
%      ROTMATRIX = FICK2ROTMATRIX(VFICK,OPTION) returns the rotation matrix
%        ROTMATRIX under the form specified by the OPTION string:
%           - 'tensor': ROTMATRIX is a 3*3*N tensor (the default)
%           - 'cell': ROTMATRIX is a 1*N cell and each cell component
%                  ROTMATRIX{1,ii} is the 3*3 rotation matrix of rotation
%                  ii
%
% Useful reference: T. Haslwanter (1995), Mathematics of three-dimensional
%    eye rotations, Vision research, 35(12), 1727-39
%
% See also ROTMATRIX2DQUAT, DQUAT2ROTMATRIX, ROTMATRIX2FICK

sv = size(vfick);
if sv == [1 3], vfick = vfick'; sv = size(vfick); end

% wrong format
if sv(1) ~= 3
    error('DualQuaternion:fick2rotMatrix:wrongsize',...
        '%d rows in the VFICK array. It should be 3.',sv(1));
end
vfick = vfick./180.*pi; % expressed in radians
n = sv(2);

% rotation matrix components (in Fick convention)
t = vfick(1,:); % Horizontal angle : + when to the left (see Haslwanter95)
f = vfick(2,:); % Vertical angle: + when down
p = vfick(3,:); % Torsional angle: + when clockwise

a11 = cos(t).*cos(f);
a22 = sin(t).*sin(f).*sin(p)+cos(t).*cos(p);
a33 = cos(f).*cos(p);
a23 = sin(t).*sin(f).*cos(p) - cos(t).*sin(p);
a32 = cos(f).*sin(p);
a31 = -sin(f);
a13 = cos(t).*sin(f).*cos(p) + sin(t).*sin(p);
a12 = cos(t).*sin(f).*sin(p)-sin(t).*cos(p);
a21 = sin(t).*cos(f);

choice = 0;
soptargin = size(varargin);
if soptargin(2) > 1
    error('DualQuaternion:fick2rotMatrix:toomanyInputs',...
        'Too Many Inputs');
elseif soptargin(2) == 1
    opt = varargin{1,1};
    if strcmp(opt,'cell') % option is 'cell'
        choice = 1;
    end
end

if choice == 0 % tensor
    rotMatrix = zeros(3,3,n);
    rotMatrix(1,1,:) = a11;
    rotMatrix(1,2,:) = a12;
    rotMatrix(1,3,:) = a13;
    rotMatrix(2,1,:) = a21;
    rotMatrix(2,2,:) = a22;
    rotMatrix(2,3,:) = a23;
    rotMatrix(3,1,:) = a31;
    rotMatrix(3,2,:) = a32;
    rotMatrix(3,3,:) = a33;
elseif choice == 1 % cell
    rotMatrix = cell(1,n);
    for ii=1:n
       rotMatrix{1,ii} = [a11(ii) a12(ii) a13(ii);...
                          a21(ii) a22(ii) a23(ii);...
                          a31(ii) a32(ii) a33(ii)];
    end
end


text:
function PARA = setMIpara(PARA,varargin)

% default parameter:
kernel  = 'cos4';
sigma   = 32;
minR    = -50;  maxR    = 300; ngvR = 8;
minT    = -50;  maxT    = 300; ngvT = 8;
entropyTol = 1e-8;

% overwrite PARA
for k=1:length(varargin)/2,
  str = sprintf('PARA=setfield(PARA,''%s'',varargin{%d});',varargin{2*k-1},2*k);
  disp(str);
  eval(str);
end;

% enpack PARA
fn = fieldnames(PARA);
for j=1:length(fn),
  str = sprintf('%s=getfield(PARA,fn{%d});',fn{j},j);
%   disp(str);
  eval(str);
end;

% check integration
tol     = 1e-2;
maxngv  = 512;

while 1,
  gvR  = linspace(minR,maxR,ngvR);
  yR   = feval(kernel,gvR,sigma);
  intR = sum(yR)*(gvR(2)-gvR(1));
  if abs(intR-1) <= tol, break; end;
  ngvR = 2*ngvR;
  if ngvR>maxngv, error('to many bins for grayvalues of R'); end;
end;
while 1,
  gvT  = linspace(minT,maxT,ngvT);
  yT   = feval(kernel,gvT,sigma);
  intT = sum(yT)*(gvT(2)-gvT(1));
  if abs(intT-1) <= tol, break; end;
  ngvT = 2*ngvT;
  if ngvT>maxngv, error('to many bins for grayvalues of T'); end;
end;

PARA.kernel = kernel;
PARA.sigma  = sigma;
PARA.minR   = minR;
PARA.maxR   = maxR;
PARA.ngvR   = ngvR;
PARA.gvR    = gvR;
PARA.minT   = minT;
PARA.maxT   = maxT;
PARA.ngvT   = ngvT;
PARA.gvT    = gvT;
PARA.entropyTol = entropyTol;

return;
% ------------------------------------------------------------------------------

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% Script for loading images LoadImages.mat and retrieve HOG

%Resizes the image to 4 times its value
%Look at the HOG features obtained to understand need for this step
%Stores features in .mat file so we dont have to run this repeatedly
function [Features] = LoadImages()
	Features = [];
	for j = 1:2
		num = num2str(j);
		str = strcat('../CIFAR10/small_data_batch_',num,'.mat');
		load(str);
		Features = vertcat(Features,data);
    end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
addpath('D:\CODE\MariusBox\SpikeDetection')

f0      = 30; % sampling frequencys

mname   = 'M150329_MP009';
datExp  = '2015-04-27';
root    = 'D:\DATA\F\';
plane_um = 0;
pix_um  = 2;
npl     = 1;

%% load the procesed data
[Ff, res] = loadPROC(root, mname, datExp, npl, plane_um, pix_um, f0);

%% ITeRATIVE OPTIMIZATION OF PARAMETERS
% load('D:\CODE\MariusBox\BigNeuralCode\results\driv1112F.mat')
load('D:\CODE\MariusBox\SpikeDetection\kernel.mat')

[dcell, Ffr, kinterp] = run_deconvolution(Ff, f0, kernel);

%%
NN = numel(dcell);
iNN = ceil(rand *NN);
clf
plot(my_conv2(Ff(:,iNN) - dcell{iNN}.B, 2, 1))
hold all
plot(conv(kinterp, Ffr(:,iNN)))
hold off

author: cortex-lab


text:
function k = prtRvUtilDiscreteKld(q,p)
% DISCRETEKLD

k = q.*(log(q)-log(p));
k(q==0) = 0;
k = sum(k(:));

end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [Simil3]=fuzzysimil3(A,B)

PA=(A(1)+2*A(2)+2*A(3)+A(4))/6;
PB=(B(1)+2*B(2)+2*B(3)+B(4))/6;
dAB=abs(PA-PB);
Simil3=1/(1+dAB);


text:
function n = exploration_noise(t, params)
n = params.noiseLevel * sin(1*t) + ...
    params.noiseLevel * sin(7*t)  + ...
    params.noiseLevel * sin(3*t);
end

author: yu-jiang

repo: radpbook

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


text:
function A_sparse = build_sparse ( )

%*****************************************************************************80
%
%% BUILD_SPARSE demonstrates how a sparse matrix can be assembled in parallel.
%
%  Discussion:
%
%    Actually, what this script demonstrates is how components of the sparse
%    matrix can be evaluated in parallel.  To create the final sparse matrix
%    data structure requires that the "sparse()" command itself be called on
%    a single processor.
%
%    Thus, this script tries to take advantage of parallel processing, assuming
%    that the evaluation of matrix entries is expensive, compared to the minor
%    cost of collecting the entries from each processor to form the final matrix.
%
%    We want to construct a rectangular sparse matrix.  This matrix can be
%    thought of as a diagonal block matrix, with each diagonal block
%    being of arbitrary size and shape.
%
%    This problem was posed by Vitor Nunes, 18 August 2011.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 August 2011
%
%  Author:
%
%    Gene Cliff
%
%  Parameters:
%
%    Output, real sparse A_sparse(*,*), the sparse matrix.
%
%  Local parameters:
%
%    Local, real A_kb(*,*), used to store one block of the matrix.
%
%    Local, integer i_blk(N_blk), the global row index just before each
%    block begins.
%
%    Local, cell array II, the global row indices for the entries
%    in each block.
%
%    Local, integer j_blk(N_blk), the global column index just before
%    each block begins.
%
%    Local, cell array JJ(N_blk,1), the global column indices for the entries
%    in each block.
%
%    Local, integer N_blk, the number of block matrices.
%
%    Local, integer row_blks(N_blk), the row dimension of each block.
%
%    Local, integer col_blks(N_blk), the column dimension of each block.
%
%    Local, cell array VV, the values of the entries in each block.
%

%
%% Set data that defines the arrangement of blocks.
%
%  The row dimension of the global matrix is the sum of row_blks,
%  and the column dimension is the sum of col_blks.  In this
%  case, our matrix has dimensions 44 by 41, and 359 nonzero entries.
%
  N_blk = 5;
  row_blks = [ 4;  7;  9; 10; 14];
  col_blks = [ 8; 10;  5; 10;  8];

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Global matrix dimension is M = %d, N = %d\n', ...
    sum ( row_blks ), sum ( col_blks ) );
%
%% Initialization.
%
  i_blk = cumsum ( row_blks(1:N_blk-1) );
  j_blk = cumsum ( col_blks(1:N_blk-1) );
  i_blk = [ 0; i_blk ];
  j_blk = [ 0; j_blk ];
%
%  The cell arrays II, JJ, and VV will contain the I, J and value
%  information for the entries in each block.
%
  II = cell ( N_blk, 1);
  JJ = cell ( N_blk, 1);
  VV = cell ( N_blk, 1);
%
%% Loop over the sub-blocks.
%
  parfor k_blk = 1 : N_blk
%
%  Here, the user would fill in the numeric values of the blocks
%  in a way that depended on the problem that was being set up.
%
%  For this demonstration, CREATE_BLOCK simply assigns arbitrary
%  values to the entries in the block.  In an actual problem, we
%  assume that the evaluation of these entries would be an expensive
%  process, so that parallel execution will result in a substantial
%  speedup.
%
    A_kb = create_block ( k_blk, row_blks(k_blk), col_blks(k_blk) );

    i_row = 1 : row_blks(k_blk);
%
%  I_BLK and J_BLK tell us how far to shift the local indices to get
%  the global row and column indices II and JJ.
%
    n_els = row_blks(k_blk) * col_blks(k_blk);

    II{k_blk}(:,1) = i_blk(k_blk) ...
      + reshape ( repmat ( i_row(:), 1, col_blks(k_blk) ), n_els, 1 );

    JJ{k_blk}(:, 1) = j_blk(k_blk) ...
      + reshape ( repmat ( 1:col_blks(k_blk), row_blks(k_blk), 1 ), n_els, 1 );

    VV{k_blk}(:, 1) = reshape ( A_kb, n_els, 1 );

  end
%
%% Now have the SPARSE command assemble the sparse matrix from the
%  data in the cell arrays.  This command does NOT execute in parallel.
%
  A_sparse = sparse ( cell2mat ( II ), cell2mat ( JJ ), cell2mat ( VV ), ...
    sum ( row_blks ), sum ( col_blks ) );

  return
end


text:
function f = olsmatrix(X,mode)

% function f = olsmatrix(X,mode)
%
% <X> is samples x parameters
% <mode> (optional) is
%   0 means normal operation
%   1 means use inv instead of \ and omit unit-length normalization.
%     the point of this mode is to reduce memory usage (i think).
%   default: 0.
%
% what we want to do is to perform OLS regression using <X>
% and obtain the parameter estimates.  this is accomplished
% by inv(X'*X)*X'*y = f*y where y is the data (samples x cases).
%
% what this function does is to return <f> which has dimensions
% parameters x samples.
%
% to ensure well-conditioning, we unit-length normalize each column
% of <X> before doing the inv.m operation.  also, we actually use
% left-slash (\), which apparently is equivalent to inv.m but faster
% and more accurate (see code for details).  if you pass <mode> as 1,
% we omit the normalization and use the inv method instead of the \ method.
%
% also, in the case that one or more regressors in <X> are all zero, then
% without special handling, then this case will result in warnings and
% NaN results.  what we do is to explicitly ensure that all-zero regressors
% are ignored and that there are all-zero rows for these regressors
% in <f>.  this makes it such that the weights estimated for these
% regressors are simply zero.
%
% see also olsmatrix2.m.
%
% history:
% 2011/06/27 - explicitly ignore 0
% 2011/03/29 - mode 1 now omits unit length normalization.
% 2011/03/28 - add <mode> input.
% 2010/08/05: now, we do not try to detect weird cases with low variance.
%             instead, we just blindly unit-length normalize each column.

% input
if ~exist('mode','var') || isempty(mode)
  mode = 0;
end

% check
assert(all(~isnan(X(:))));

% deal with degenerate regressors
good = ~all(X==0,1);

% initialize result
f = zeros(size(X,2),size(X,1));

% do it
switch mode
case 0
  [X,len] = unitlength(X(:,good),1,[],0);
  temp = diag(1./len)*((X'*X)\X');  % hm, suggested by M-Lint
case 1
  temp = inv(X(:,good)'*X(:,good))*X(:,good)';
end

% return
if any(good)
  f(good,:) = temp;
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab