text:
%kksampline 'Sample a data object along an arbitrary line'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros ksampline.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'input data object'
% Integer: n 'Number of samples (including endpoints)', default: 1000: 'number of samples (including endpoints)'
% Integer: nc 'Minimum required neighbor count', default: 0: 'minimum required neighbor count'
% OutputFile: o 'Output', required: 'output data object'
%
% Example: o = kksampline(i, {'i','';'n',1000;'nc',0;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% ksampline - Sample a data object along an arbitrary line
%
%  DESCRIPTION
% .I ksampline
% is used to sample the data in an object along an arbitrary line and at
% arbitrary intervals. The value at each sample point is obtained by
% an inverse distance weighting scheme.
%
% The first sample point is precisely at the starting coordinate and
% the final sample is precisely at the ending coordinate.
%
% The line may traverse any portion of the 5D data space. If a particular
% sample on the line lies outside the data space then the values will be
% interpolated using a padding value of zero.
%
% The output object will be of the same data type as the input object. If
% a map is present in the input object, then the value data is sent through the
% map before processing. The output data will thus have no map.
%
% Mask data, if present, is used to help control the interpolation.
% Value points with a corresponding mask value of zero will not be included
% in the interpolated value. However, a minimum number of neighboring
% values with a valid mask are required in order to produce a valid output
% value; non-valid output values are marked by a zero mask. If the mask
% segment is not present, then all value points in the source object are assumed
% to be valid.
%
% The actual minimum
% required neighbor count is determined by the setting of the "nc" argument.
% If nc is zero, then the minimum required neighbor count will be computed from
% the dimensionality of the source data set. For nc=0, the minimum required
% neighbor count is 0.5*(2^dim)+1 where dim is the dimensionality of
% the source object (1 for a line, 2 for a plane, 3 for a volume, etc).
% If nc is non-zero, then the minimum required neighbor count is simply nc.
% Beware of non-intuitive results when interpolating across mask wise non-convex
% parts of the data object. The default minimum required neighbor count
% gives "reasonable" answers in these cases. Reducing the minimum required
% neighbor count can provide answers that are less and less meaningful.
%
% Note that ksampline samples the data as if it were using a delta function.
% This implies that you may get significantly different
% results depending on the position of each sample location. You may wish to
% lowpass the source data before sampling with ksampline to avoid aliasing.
% This can be particularly important when dealing with "holey data" or when
% the sampling involves degenerate geometry (such as sampling a plane object
% with a sampling line that does not lie in the plane).
%
% ksampline is implemented using "point\fR data access. A point data access
% is done for each sample in the line. This implies that a large number of samples
% may require quite some time to compute, particularly if operating on a
% higher dimension data set.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kksampline(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kksampline(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'n', 1000;'nc', 0;'o', '__output'};
maxval={0,0,32,0};
minval={0,0,0,0};
istoggle=[0,1,1,0];
was_set=istoggle * 0;
paramtype={'InputFile','Integer','Integer','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'ksampline"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function varargout=overlapadd(varargin)
%
% For calling details please see v_overlapadd.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_overlapadd.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: overlapadd.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_overlapadd(varargin{:});
else
    v_overlapadd(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [z,dz,ymu,ys,fmu,fs,fpi] = acqNegPI(xi,target,gpstruct,optimState,grad_flag)
%ACQNEGPI Acquisition function for (negative) probability of improvement.

if nargin < 5 || isempty(grad_flag); grad_flag = false; end

n = size(xi,1);
Nhyp = numel(gpstruct.hyp);

if grad_flag && n > 1
    error('acqNegPI:gradient', ...
        'Gradient of acquisition function is provided only at one test point XI (row vector).');
end

if grad_flag
    [ymu,ys2,fmu,fs2,hypw,dymu,dys2,dfmu,dfs2] = gppred(xi,gpstruct,'central');
else
    [ymu,ys2,fmu,fs2,hypw] = gppred(xi,gpstruct);
end
fs = sqrt(fs2);
ys = sqrt(ys2);

% Negative probability of improvement
gammaz = (target - fmu)./fs;
z = -0.5*erfc(-gammaz/sqrt(2));

try
    z = sum(bsxfun(@times,hypw(~isnan(hypw)),z(~isnan(hypw),:)),1);
catch
    z = Inf(1,n);
    dz = NaN(n,size(xi,2));
    return;
end

if grad_flag
    % Gradient of probability of improvement
    dfs = 0.5*dfs2./fs;
    dgammaz = -(dfmu.*fs + (target - fmu).*dfs)./fs2;
    dz = 0.5*dgammaz/sqrt(2)*(-2*exp(-gammaz.^2/2)/sqrt(pi));

    dz = sum(bsxfun(@times,hypw,dz(~isnan(hypw),:)),1);
else
    dz = NaN(n,size(xi,2));     % Gradient not estimated
end

end

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
function [vn,ins_align_flag]=tdcp2vel(rtk,nav,cur_obs,old_obs)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2020-2025, by Kai Chen, All rights reserved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global glc
vn=zeros(1,3); ins_align_flag=0;
dt=timediff(cur_obs(1).time,old_obs(1).time);

cur_sv=satposs(cur_obs,nav,rtk.opt.sateph);
old_sv=satposs(old_obs,nav,rtk.opt.sateph);

cur_rr=rtk.sol.pos;  [cur_pos,Cne]=xyz2blh(cur_rr);
old_rr=rtk.sol.pos;  old_pos=xyz2blh(old_rr);

n_cur=size(cur_obs,1);
n_old=size(old_obs,1);
v=zeros(n_cur,1); H=zeros(n_cur,4); P=zeros(n_cur,n_cur);

nv=0;
for i=1:n_cur

    cur_obsi=cur_obs(i);cur_svi=cur_sv(i);
    lam=nav.lam(cur_obsi.sat,:);
    [sys,~]=satsys(cur_obsi.sat);
    if rtk.mask(sys)==0,continue;end

    if lam(1)==0,continue;end

    old_idx=0;
    for j=1:n_old
        if cur_obsi.sat==old_obs(j).sat
            old_idx=j;break;
        end
    end
    if old_idx==0,continue;end

    old_obsi=old_obs(old_idx); old_svi=old_sv(old_idx);

    if cur_svi.svh~=0||old_svi.svh~=0||norm(cur_svi.pos)<=0||norm(old_svi.pos)<=0
        continue;
    end
    if cur_obsi.L(1)==0||old_obsi.L(1)==0,continue;end

    [cur_range,cur_LOS]=geodist(cur_svi.pos,cur_rr');
    cur_azel=satazel(cur_pos,cur_LOS);
    [old_range,old_LOS]=geodist(old_svi.pos,old_rr');
    old_azel=satazel(old_pos',old_LOS); %#ok

    tdcp_obs=cur_obsi.L(1)*lam(1)-old_obsi.L(1)*lam(1);
    tdcp_pre=cur_range-old_range;

    v(nv+1,1)=tdcp_obs-tdcp_pre+glc.CLIGHT*(cur_svi.dts-old_svi.dts);

    H(nv+1,:)=[-cur_LOS,1];

    P(nv+1,nv+1)=(2*1^2/sin(cur_azel(2))^2)^-1;

    nv=nv+1;

end

if nv<4,return;end
if nv<n_cur
    v(nv+1:end,:)=[]; H(nv+1:end,:)=[]; P(nv+1:end,:)=[];P(:,nv+1:end)=[];
end

%robust
nv=size(v,1);exc=zeros(nv,1);
v0=zeros(nv,1); H0=zeros(nv,4); P0=zeros(nv,nv); n=0;
for i=1:nv
    vi=abs(v(i));
    v_other=abs(v);v_other(i)=[];
    ave_v=sum(abs(v_other))/(nv-1);
    % the average distance from the population value
    ave_distance=sum(abs(vi-v_other))/(nv-1);
    if ave_distance>3*ave_v||(ave_distance>100&&ave_distance>1.2*ave_v)
        exc(i)=1;
    end
end
for i=1:nv
    if exc(i)==1,continue;end
    v0(n+1,1)=v(i);
    H0(n+1,:)=H(i,:);
    P0(n+1,n+1)=P(i,i);
    n=n+1;
end
if size(v0,1)<4,return;end
if n<nv
    v0(n+1:end,:)=[]; H0(n+1:end,:)=[]; P0(n+1:end,:)=[];P0(:,n+1:end)=[];
end

[dx,~]=least_square(v0,H0,P0);

ve=dx(1:3)/dt;
vn=(Cne*ve)';

if dot(vn,vn)>3
    ins_align_flag=1;
end

return

repo: GINav


text:
classdef tblock < matlab.unittest.TestCase
    % tblock   Unit tests for transformer.layer.block

    % Copyright 2020 The MathWorks, Inc.

    properties(Constant, Access=private)
        block = @gpt2.layer.block
    end

    properties(TestParameter)
        Input = struct(...
            'Scalar', 1,...
            'Vector', 1:5,...
            'Matrix', reshape(1:6,[3,2]))
    end

    methods(Test)
        function outputHasInputSize(test,Input)
            % The block is simply a composition of other layers. Simply
            % verify the output of a block is the same size as the input,
            % and as such the blocks can be stacked.
            x = dlarray(Input);
            C = size(Input,1);
            weights = test.randomWeights(C);
            hyperParameters.NumHeads = 1;
            y = test.block(x,[],weights,hyperParameters);
            test.verifySize(y,size(x));
        end

        function outputHasInputSizeWithPasts(test,Input)
            % As above but using "pasts" - a concatenation of key and value
            % matrices.
            x = dlarray(Input);
            C = size(Input,1);
            weights = test.randomWeights(C);
            hyperParameters.NumHeads = 1;
            % Provide a fake past of sequence length 1
            K_fake = dlarray(rand(C,1));
            V_fake = dlarray(rand(C,1));
            past = cat(5,K_fake,V_fake);
            [y,present] = test.block(x,past,weights,hyperParameters);
            test.verifySize(y,size(x));
            % The size of presents is the size of past except the sequence
            % dimension gets extended by the sequence length of y
            exp_present_size = size(past);
            exp_present_size(2) = exp_present_size(2)+size(y,2);
            test.verifySize(present,exp_present_size);
        end
    end

    methods(Access=private)
        function weights = randomWeights(test,C)
            % C is num features, or latent dimension of the block
            g1 = dlarray(rand(C,1));
            b1 = dlarray(rand(C,1));
            g2 = dlarray(rand(C,1));
            b2 = dlarray(rand(C,1));
            W_A1 = dlarray(rand(3*C,C));
            W_A2 = dlarray(rand(C));
            b_A1 = dlarray(rand(3*C,1));
            b_A2 = dlarray(rand(C,1));
            W_P1 = dlarray(rand(C));
            b_P1 = dlarray(rand(C,1));
            W_P2 = dlarray(rand(C));
            b_P2 = dlarray(rand(C,1));
            weights = test.prepareBlockWeightsStruct(g1,b1,W_A1,b_A1,W_A2,b_A2,g2,b2,W_P1,b_P1,W_P2,b_P2);
        end

        function s = prepareBlockWeightsStruct(test,g1,b1,W_A1,b_A1,W_A2,b_A2,g2,b2,W_P1,b_P1,W_P2,b_P2)
            % Merge various structs that have the appropriate weight naming
            % syntax.
            s_ln = test.prepareLayerNormWeightsStruct(g1,b1,g2,b2);
            s_attn = test.prepareAttentionWeightsStruct(W_A1,b_A1,W_A2,b_A2);
            s_mlp = test.prepareMLPWeightsStruct(W_P1,b_P1,W_P2,b_P2);
            c = {s_ln,s_attn,s_mlp};
            fn = cellfun(@fieldnames,c,'UniformOutput',false);
            fn = cat(1,fn{:});
            fv = cellfun(@struct2cell,c,'UniformOutput',false);
            fv = cat(1,fv{:});
            s = struct();
            for i = 1:numel(fn)
                s.(fn{i}) = fv{i};
            end
        end

        function s = prepareAttentionWeightsStruct(~,W1,b1,W2,b2)
            % Prepare a struct compatible with the weights input of
            % attention. These are for the fully connected layers.
            s = struct(...
                'attn_c_attn_w_0',W1,...
                'attn_c_attn_b_0',b1,...
                'attn_c_proj_w_0',W2,...
                'attn_c_proj_b_0',b2);
        end

        function s = prepareLayerNormWeightsStruct(~,g1,b1,g2,b2)
            % Prepare a struct of weights compatible with the two layer
            % norm calls in block
            s = struct(...
                'ln_1_g_0',g1,...
                'ln_1_b_0',b1,...
                'ln_2_g_0',g2,...
                'ln_2_b_0',b2);
        end

        function s = prepareMLPWeightsStruct(~,W1,b1,W2,b2)
            % Create a struct of weights to be consumed by
            % transformer.layer.multiLayerPerceptron
            s = struct(...
                'mlp_c_fc_w_0',W1,...
                'mlp_c_fc_b_0',b1,...
                'mlp_c_proj_w_0',W2,...
                'mlp_c_proj_b_0',b2);
        end
    end
end

author: matlab-deep-learning

repo: transformer-models


text:
%%  kmeans algorithm for an image
%---input---------------------------------------------------------
%   Y: 3D image
%   k: number of clusters
%   g: number of GMM components
%---output--------------------------------------------------------
%   X: 3D labels
%   GMM: Gaussian mixture model parameters

function [X GMM]=image_kmeans(Y,k,g)
[m n l]=size(Y);
y=Y(:);
x=kmeans(y,k);
X=reshape(x,[m n l]);

GMM=get_GMM(X,Y,g);

author: Sable

repo: mcbench-benchmarks


text:
function [alpha,b,eta,zeta,kappa] = solve_hkl_logistic_dual_fast(eta,data,Y,lambda,varargin);

% solve the milasso for a given regularization parameter
n = size(data.X,1);
p =length(data.affinity);
pX = size(data.X,2);
affs = data.affinity;
weights = data.weights;

% optional parameters
mingap = 1e-3;
display = 0;
gapeta = 1e-3;

% READ OPTIONAL PARAMETERS
args = varargin;
nargs = length(args);
for i=1:2:nargs
	switch args{i},
		case 'mingap',        mingap = args{i+1};
		case 'gapeta',        gapeta = args{i+1};
		case 'display',        display = args{i+1};
	end
end

if isempty(eta)
	eta = 1/p ./ weights.^2;
end
x = eta .* weights.^2;
x = max(x,0);
x = x / sum(x);

if display==1
	display = Inf;
end
beta_ARMIJO = .25 ;
sigma_ARMIJO = .1 ;
kmax_ARMIJO = 20;
kmax = 200;
tol_dx = 1e-12;
k = 1;
alpha=1;

% % optimization parameters for the dual problem
% optparam.tol = 1e-16;
% optparam.kmax = 200;
% optparam.tol_df = 1e-20;
% optparam.display= 0;
fx0 = Inf;

% solve the regular kernel learning problem
X = data.X;	% no need to center here for logistic
mY = mean(Y);
G = zeros(size(data.X,2),p);
for i=1:p
	G(data.groups{i},i) = 1;
end

% optimization parameter for logistic regression
optparam.tol = 1e-16;
optparam.kmax = 50;
optparam.tol_df = 1e-14;
optparam.display= 0;

%nu =  1e-6; % added ridge on kernels
while k<=kmax;
	x = max(x,0);
	x = x / sum(x);

	% compute value of function
	eta = ( x*(1-gapeta) + gapeta/p )./ weights.^2 ;
	zeta = zeros(p,1);
	for i=1:p
		zeta(affs{i}) = zeta(affs{i}) + 1./eta(i);
	end
	zeta = 1./zeta;
		reweight = G * sqrt(zeta);

	% solve the regular kernel learning problem in the primal
	if k==1
		% first iteration -> initialize with square loss
		% otherwise, take previous iteration weights
		temp = ( reweight .* ( X' * (Y - mY ) ) );
		wt = (  (reweight*reweight') .* ( X' * X ) + n * lambda * eye(pX) ) \   temp ;
	end

	Xr = repmat(reweight',n,1) .* X;
	[wt,exit_parameters] = minimize_newton(wt,@logistic_cost,optparam, Xr,Y,lambda);
	fx = logistic_cost(wt,Xr,Y,lambda);

	gradient_zeta = zeros(p,1);
	for i=1:p
		gradient_zeta(i) = - lambda/2 *  norm(wt(data.groups{i}))^2 / zeta(i);
	end

	gradient_eta = zeros(p,1);
	for i=1:p;
		gradient_eta(i) = sum( gradient_zeta(affs{i}) .* zeta(affs{i}).^2 ) / eta(i)^2;
	end
	gradient = (1-gapeta) * gradient_eta ./ weights.^2;

 alpha = 1;
	% check optimality condition
	optcond = max(x .* ( gradient - min(gradient)));
	if max(x .* ( gradient - min(gradient))) < mingap, break; end

	fx0=fx;

	% projected gradient
	xbar = x -   gradient;
	d = - gradient;
	x0=x;
	ialpha = 1;
	while ialpha<kmax_ARMIJO
		xalpha = simplex_projection(x + d * alpha);
 		% fxalpha = feval(fun,xalpha,varargin{:});
		eta = ( xalpha*(1-gapeta) + gapeta/p )./ weights.^2 ;
		zeta = zeros(p,1);
		for i=1:p
			zeta(affs{i}) = zeta(affs{i}) + 1./eta(i);
		end
		zeta = 1./zeta;

		reweight = G * sqrt(zeta);
		Xr = repmat(reweight',n,1) .* X;
		[wt0,exit_parameters] = minimize_newton(wt,@logistic_cost,optparam, Xr,Y,lambda);
		fxalpha = logistic_cost(wt0,Xr,Y,lambda);

 		if fx - fxalpha >= - sigma_ARMIJO * gradient' * ( xalpha - x );
 			% start to go up
			while ( fx - fxalpha >= - sigma_ARMIJO * gradient' * ( xalpha - x ) ) & ialpha<kmax_ARMIJO;
				alpha = alpha / beta_ARMIJO;
 				xalpha = simplex_projection(x + d * alpha);
				%fxalpha = feval(fun,xalpha,varargin{:});
				eta = ( xalpha*(1-gapeta) + gapeta/p )./ weights.^2 ;
				zeta = zeros(p,1);
				for i=1:p
					zeta(affs{i}) = zeta(affs{i}) + 1./eta(i);
				end
				zeta = 1./zeta;

				reweight = G * sqrt(zeta);
				Xr = repmat(reweight',n,1) .* X;
				[wt0,exit_parameters] = minimize_newton(wt,@logistic_cost,optparam, Xr,Y,lambda);
				fxalpha = logistic_cost(wt0,Xr,Y,lambda);

				% fprintf('+');
				ialpha = ialpha+1;
			end
			alpha = alpha * beta_ARMIJO;
			xalpha = simplex_projection(x + d * alpha);
			x = xalpha;
			break;
		else
			alpha = alpha * beta_ARMIJO;
 		end
		ialpha = ialpha+1;
	end
 	if ialpha==kmax_ARMIJO, alpha=1; end

	if ~isinf(display)
		if display==1,
			fprintf('k=%d - f=%f - armijo=%d - dx=%e - gap=%f\n',k,fx,ialpha,norm(x-x0),optcond);
		end
		if display>1,
			if mod(k,display)==1
				fprintf('k=%d - f=%f - armijo=%d - dx=%e - gap=%f\n',k,fx,ialpha,norm(x-x0),optcond);

			end
		end
	end
	if norm(x-x0)<tol_dx,    break; end

	if sum(x) > 1+1e-1, keyboard; end

	k=k+1;
end

% if k==kmax+1
%     keyboard;
% end

eta = ( x*(1-gapeta) + gapeta/p )./ weights.^2 ;
zeta = zeros(p,1);
for i=1:p
	zeta(affs{i}) = zeta(affs{i}) + 1./eta(i);
end
zeta = 1./zeta;
reweight = G * sqrt(zeta);
Xr = repmat(reweight',n,1) .* X;
[wt,exit_parameters] = minimize_newton(wt,@logistic_cost,optparam, Xr,Y,lambda);

pred = data.X * ( reweight.* wt);
b= center_logistic(pred,Y);

alpha = 1/(n*lambda) * ( Y - 1./( 1 + exp( - pred - b ) ) );

kappa = zeros(p,p);
for i=1:p
	kappa(i,affs{i}) = zeta(affs{i})./eta(i);
end

if isinf(display)
	fprintf('k=%d - f=%f -  gap=%f\n',k,fx,optcond);
end
% if optcond> 1,
% keyboard
% end

repo: REST


text:
% The COBRAToolbox: testFASTCC.m
%
% Purpose:
%     - test FASTCC algorithm
%
% Authors:
%     - Original file: Thomas Pfau, May 2016
%     - Fix by @fplanes July 2017
%     - CI integration: Laurent Heirendt July 2017
%

% save the current path
currentDir = pwd;

% initialize the test
fileDir = fileparts(which('testFASTCC'));
cd(fileDir);

%load a model
model = getDistributedModel('ecoli_core_model.mat');
% create a model with fructose
modelWithFru = changeRxnBounds(model,'EX_fru(e)',-100,'l');

% set paarmeters
epsilon = 1e-4;
printLevel = 2;
modeFlag = 0;

%by default, fructose updatek, fumarate uptake and corresponding reactions
%cannot be used.
inactives = {'EX_fru(e)', 'EX_fum(e)','EX_gln_L(e)', 'EX_mal_L(e)', 'FRUpts2', 'FUMt2_2', 'GLNabc','MALt2_2' };
fructoseRelated = {'EX_fru(e)', 'FRUpts2'};

% define the solver packages to be used to run this test
solverPkgs = prepareTest('needsLP',true);

for k = 1:length(solverPkgs.LP)
    % change the COBRA solver (LP)
    solverOK = changeCobraSolver(solverPkgs.LP{k}, 'LP', 0);
    fprintf('   Testing FASTCC using %s ... \n', solverPkgs.LP{k});
    A = fastcc(model, epsilon, printLevel,modeFlag);
    assert(isempty(setxor(setdiff(model.rxns,model.rxns(A)),inactives)));
    % Open up the Fructose channel
    A = fastcc(modelWithFru, epsilon, printLevel,modeFlag);
    % now, everything from the inactives except the fructose reactions are
    % not in A
    assert(isempty(setxor(setdiff(modelWithFru.rxns,modelWithFru.rxns(A)),setdiff(inactives,fructoseRelated))));
    % output a success message
    fprintf('Done.\n');
end

% change the directory
cd(currentDir)

author: opencobra

repo: cobratoolbox


text:
function varargout = colorspace(Conversion,varargin)
%COLORSPACE  Convert a color image between color representations.
%   B = COLORSPACE(S,A) converts the color representation of image A
%   where S is a string specifying the conversion.  S tells the
%   source and destination color spaces, S = 'dest<-src', or
%   alternatively, S = 'src->dest'.  Supported color spaces are
%
%     'RGB'              R'G'B' Red Green Blue (ITU-R BT.709 gamma-corrected)
%     'YPbPr'            Luma (ITU-R BT.601) + Chroma
%     'YCbCr'/'YCC'      Luma + Chroma ("digitized" version of Y'PbPr)
%     'YUV'              NTSC PAL Y'UV Luma + Chroma
%     'YIQ'              NTSC Y'IQ Luma + Chroma
%     'YDbDr'            SECAM Y'DbDr Luma + Chroma
%     'JPEGYCbCr'        JPEG-Y'CbCr Luma + Chroma
%     'HSV'/'HSB'        Hue Saturation Value/Brightness
%     'HSL'/'HLS'/'HSI'  Hue Saturation Luminance/Intensity
%     'XYZ'              CIE XYZ
%     'Lab'              CIE L*a*b* (CIELAB)
%     'Luv'              CIE L*u*v* (CIELUV)
%     'Lch'              CIE L*ch (CIELCH)
%
%  All conversions assume 2 degree observer and D65 illuminant.  Color
%  space names are case insensitive.  When R'G'B' is the source or
%  destination, it can be omitted. For example 'yuv<-' is short for
%  'yuv<-rgb'.
%
%  MATLAB uses two standard data formats for R'G'B': double data with
%  intensities in the range 0 to 1, and uint8 data with integer-valued
%  intensities from 0 to 255.  As MATLAB's native datatype, double data is
%  the natural choice, and the R'G'B' format used by colorspace.  However,
%  for memory and computational performance, some functions also operate
%  with uint8 R'G'B'.  Given uint8 R'G'B' color data, colorspace will
%  first cast it to double R'G'B' before processing.
%
%  If A is an Mx3 array, like a colormap, B will also have size Mx3.
%
%  [B1,B2,B3] = COLORSPACE(S,A) specifies separate output channels.
%  COLORSPACE(S,A1,A2,A3) specifies separate input channels.

% Pascal Getreuer 2005-2006

%%% Input parsing %%%
if nargin < 2, error('Not enough input arguments.'); end
[SrcSpace,DestSpace] = parse(Conversion);

if nargin == 2
   Image = varargin{1};
elseif nargin >= 3
   Image = cat(3,varargin{:});
else
   error('Invalid number of input arguments.');
end

FlipDims = (size(Image,3) == 1);

if FlipDims, Image = permute(Image,[1,3,2]); end
if ~isa(Image,'double'), Image = double(Image)/255; end
if size(Image,3) ~= 3, error('Invalid input size.'); end

SrcT = gettransform(SrcSpace);
DestT = gettransform(DestSpace);

if ~ischar(SrcT) & ~ischar(DestT)
   % Both source and destination transforms are affine, so they
   % can be composed into one affine operation
   T = [DestT(:,1:3)*SrcT(:,1:3),DestT(:,1:3)*SrcT(:,4)+DestT(:,4)];
   Temp = zeros(size(Image));
   Temp(:,:,1) = T(1)*Image(:,:,1) + T(4)*Image(:,:,2) + T(7)*Image(:,:,3) + T(10);
   Temp(:,:,2) = T(2)*Image(:,:,1) + T(5)*Image(:,:,2) + T(8)*Image(:,:,3) + T(11);
   Temp(:,:,3) = T(3)*Image(:,:,1) + T(6)*Image(:,:,2) + T(9)*Image(:,:,3) + T(12);
   Image = Temp;
elseif ~ischar(DestT)
   Image = rgb(Image,SrcSpace);
   Temp = zeros(size(Image));
   Temp(:,:,1) = DestT(1)*Image(:,:,1) + DestT(4)*Image(:,:,2) + DestT(7)*Image(:,:,3) + DestT(10);
   Temp(:,:,2) = DestT(2)*Image(:,:,1) + DestT(5)*Image(:,:,2) + DestT(8)*Image(:,:,3) + DestT(11);
   Temp(:,:,3) = DestT(3)*Image(:,:,1) + DestT(6)*Image(:,:,2) + DestT(9)*Image(:,:,3) + DestT(12);
   Image = Temp;
else
   Image = feval(DestT,Image,SrcSpace);
end

%%% Output format %%%
if nargout > 1
   varargout = {Image(:,:,1),Image(:,:,2),Image(:,:,3)};
else
   if FlipDims, Image = permute(Image,[1,3,2]); end
   varargout = {Image};
end

return;

function [SrcSpace,DestSpace] = parse(Str)
% Parse conversion argument

if isstr(Str)
   Str = lower(strrep(strrep(Str,'-',''),' ',''));
   k = find(Str == '>');

   if length(k) == 1         % Interpret the form 'src->dest'
      SrcSpace = Str(1:k-1);
      DestSpace = Str(k+1:end);
   else
      k = find(Str == '<');

      if length(k) == 1      % Interpret the form 'dest<-src'
         DestSpace = Str(1:k-1);
         SrcSpace = Str(k+1:end);
      else
         error(['Invalid conversion, ''',Str,'''.']);
      end
   end

   SrcSpace = alias(SrcSpace);
   DestSpace = alias(DestSpace);
else
   SrcSpace = 1;             % No source pre-transform
   DestSpace = Conversion;
   if any(size(Conversion) ~= 3), error('Transformation matrix must be 3x3.'); end
end
return;

function Space = alias(Space)
Space = strrep(Space,'cie','');

if isempty(Space)
   Space = 'rgb';
end

switch Space
case {'ycbcr','ycc'}
   Space = 'ycbcr';
case {'hsv','hsb'}
   Space = 'hsv';
case {'hsl','hsi','hls'}
   Space = 'hsl';
case {'rgb','yuv','yiq','ydbdr','ycbcr','jpegycbcr','xyz','lab','luv','lch'}
   return;
end
return;

function T = gettransform(Space)
% Get a colorspace transform: either a matrix describing an affine transform,
% or a string referring to a conversion subroutine
switch Space
case 'ypbpr'
   T = [0.299,0.587,0.114,0;-0.1687367,-0.331264,0.5,0;0.5,-0.418688,-0.081312,0];
case 'yuv'
   % R'G'B' to NTSC/PAL YUV
   % Wikipedia: http://en.wikipedia.org/wiki/YUV
   T = [0.299,0.587,0.114,0;-0.147,-0.289,0.436,0;0.615,-0.515,-0.100,0];
case 'ydbdr'
   % R'G'B' to SECAM YDbDr
   % Wikipedia: http://en.wikipedia.org/wiki/YDbDr
   T = [0.299,0.587,0.114,0;-0.450,-0.883,1.333,0;-1.333,1.116,0.217,0];
case 'yiq'
   % R'G'B' in [0,1] to NTSC YIQ in [0,1];[-0.595716,0.595716];[-0.522591,0.522591];
   % Wikipedia: http://en.wikipedia.org/wiki/YIQ
   T = [0.299,0.587,0.114,0;0.595716,-0.274453,-0.321263,0;0.211456,-0.522591,0.311135,0];
case 'ycbcr'
   % R'G'B' (range [0,1]) to ITU-R BRT.601 (CCIR 601) Y'CbCr
   % Wikipedia: http://en.wikipedia.org/wiki/YCbCr
   % Poynton, Equation 3, scaling of R'G'B to Y'PbPr conversion
   T = [65.481,128.553,24.966,16;-37.797,-74.203,112.0,128;112.0,-93.786,-18.214,128];
case 'jpegycbcr'
   % Wikipedia: http://en.wikipedia.org/wiki/YCbCr
   T = [0.299,0.587,0.114,0;-0.168736,-0.331264,0.5,0.5;0.5,-0.418688,-0.081312,0.5]*255;
case {'rgb','xyz','hsv','hsl','lab','luv','lch'}
   T = Space;
otherwise
   error(['Unknown color space, ''',Space,'''.']);
end
return;

function Image = rgb(Image,SrcSpace)
% Convert to Rec. 709 R'G'B' from 'SrcSpace'
switch SrcSpace
case 'rgb'
   return;
case 'hsv'
   % Convert HSV to R'G'B'
   Image = huetorgb((1 - Image(:,:,2)).*Image(:,:,3),Image(:,:,3),Image(:,:,1));
case 'hsl'
   % Convert HSL to R'G'B'
   L = Image(:,:,3);
   Delta = Image(:,:,2).*min(L,1-L);
   Image = huetorgb(L-Delta,L+Delta,Image(:,:,1));
case {'xyz','lab','luv','lch'}
   % Convert to CIE XYZ
   Image = xyz(Image,SrcSpace);
   % Convert XYZ to RGB
   T = [3.240479,-1.53715,-0.498535;-0.969256,1.875992,0.041556;0.055648,-0.204043,1.057311];
   R = T(1)*Image(:,:,1) + T(4)*Image(:,:,2) + T(7)*Image(:,:,3);  % R
   G = T(2)*Image(:,:,1) + T(5)*Image(:,:,2) + T(8)*Image(:,:,3);  % G
   B = T(3)*Image(:,:,1) + T(6)*Image(:,:,2) + T(9)*Image(:,:,3);  % B
   % Desaturate and rescale to constrain resulting RGB values to [0,1]
   AddWhite = -min(min(min(R,G),B),0);
   Scale = max(max(max(R,G),B)+AddWhite,1);
   R = (R + AddWhite)./Scale;
   G = (G + AddWhite)./Scale;
   B = (B + AddWhite)./Scale;
   % Apply gamma correction to convert RGB to Rec. 709 R'G'B'
   Image(:,:,1) = gammacorrection(R);  % R'
   Image(:,:,2) = gammacorrection(G);  % G'
   Image(:,:,3) = gammacorrection(B);  % B'
otherwise  % Conversion is through an affine transform
   T = gettransform(SrcSpace);
   temp = inv(T(:,1:3));
   T = [temp,-temp*T(:,4)];
   R = T(1)*Image(:,:,1) + T(4)*Image(:,:,2) + T(7)*Image(:,:,3) + T(10);
   G = T(2)*Image(:,:,1) + T(5)*Image(:,:,2) + T(8)*Image(:,:,3) + T(11);
   B = T(3)*Image(:,:,1) + T(6)*Image(:,:,2) + T(9)*Image(:,:,3) + T(12);
   AddWhite = -min(min(min(R,G),B),0);
   Scale = max(max(max(R,G),B)+AddWhite,1);
   R = (R + AddWhite)./Scale;
   G = (G + AddWhite)./Scale;
   B = (B + AddWhite)./Scale;
   Image(:,:,1) = R;
   Image(:,:,2) = G;
   Image(:,:,3) = B;
end

% Clip to [0,1]
Image = min(max(Image,0),1);
return;

function Image = xyz(Image,SrcSpace)
% Convert to CIE XYZ from 'SrcSpace'
WhitePoint = [0.950456,1,1.088754];

switch SrcSpace
case 'xyz'
   return;
case 'luv'
   % Convert CIE L*uv to XYZ
   WhitePointU = (4*WhitePoint(1))./(WhitePoint(1) + 15*WhitePoint(2) + 3*WhitePoint(3));
   WhitePointV = (9*WhitePoint(2))./(WhitePoint(1) + 15*WhitePoint(2) + 3*WhitePoint(3));
   L = Image(:,:,1);
   Y = (L + 16)/116;
   Y = invf(Y)*WhitePoint(2);
   U = Image(:,:,2)./(13*L + 1e-6*(L==0)) + WhitePointU;
   V = Image(:,:,3)./(13*L + 1e-6*(L==0)) + WhitePointV;
   Image(:,:,1) = -(9*Y.*U)./((U-4).*V - U.*V);                  % X
   Image(:,:,2) = Y;                                             % Y
   Image(:,:,3) = (9*Y - (15*V.*Y) - (V.*Image(:,:,1)))./(3*V);  % Z
case {'lab','lch'}
   Image = lab(Image,SrcSpace);
   % Convert CIE L*ab to XYZ
   fY = (Image(:,:,1) + 16)/116;
   fX = fY + Image(:,:,2)/500;
   fZ = fY - Image(:,:,3)/200;
   Image(:,:,1) = WhitePoint(1)*invf(fX);  % X
   Image(:,:,2) = WhitePoint(2)*invf(fY);  % Y
   Image(:,:,3) = WhitePoint(3)*invf(fZ);  % Z
otherwise   % Convert from some gamma-corrected space
   % Convert to Rec. 701 R'G'B'
   Image = rgb(Image,SrcSpace);
   % Undo gamma correction
   R = invgammacorrection(Image(:,:,1));
   G = invgammacorrection(Image(:,:,2));
   B = invgammacorrection(Image(:,:,3));
   % Convert RGB to XYZ
   T = inv([3.240479,-1.53715,-0.498535;-0.969256,1.875992,0.041556;0.055648,-0.204043,1.057311]);
   Image(:,:,1) = T(1)*R + T(4)*G + T(7)*B;  % X
   Image(:,:,2) = T(2)*R + T(5)*G + T(8)*B;  % Y
   Image(:,:,3) = T(3)*R + T(6)*G + T(9)*B;  % Z
end
return;

function Image = hsv(Image,SrcSpace)
% Convert to HSV
Image = rgb(Image,SrcSpace);
V = max(Image,[],3);
S = (V - min(Image,[],3))./(V + (V == 0));
Image(:,:,1) = rgbtohue(Image);
Image(:,:,2) = S;
Image(:,:,3) = V;
return;

function Image = hsl(Image,SrcSpace)
% Convert to HSL
switch SrcSpace
case 'hsv'
   % Convert HSV to HSL
   MaxVal = Image(:,:,3);
   MinVal = (1 - Image(:,:,2)).*MaxVal;
   L = 0.5*(MaxVal + MinVal);
   temp = min(L,1-L);
   Image(:,:,2) = 0.5*(MaxVal - MinVal)./(temp + (temp == 0));
   Image(:,:,3) = L;
otherwise
   Image = rgb(Image,SrcSpace);  % Convert to Rec. 701 R'G'B'
   % Convert R'G'B' to HSL
   MinVal = min(Image,[],3);
   MaxVal = max(Image,[],3);
   L = 0.5*(MaxVal + MinVal);
   temp = min(L,1-L);
   S = 0.5*(MaxVal - MinVal)./(temp + (temp == 0));
   Image(:,:,1) = rgbtohue(Image);
   Image(:,:,2) = S;
   Image(:,:,3) = L;
end
return;

function Image = lab(Image,SrcSpace)
% Convert to CIE L*a*b* (CIELAB)
WhitePoint = [0.950456,1,1.088754];

switch SrcSpace
case 'lab'
   return;
case 'lch'
   % Convert CIE L*CH to CIE L*ab
   C = Image(:,:,2);
   Image(:,:,2) = cos(Image(:,:,3)*pi/180).*C;  % a*
   Image(:,:,3) = sin(Image(:,:,3)*pi/180).*C;  % b*
otherwise
   Image = xyz(Image,SrcSpace);  % Convert to XYZ
   % Convert XYZ to CIE L*a*b*
   X = Image(:,:,1)/WhitePoint(1);
   Y = Image(:,:,2)/WhitePoint(2);
   Z = Image(:,:,3)/WhitePoint(3);
   fX = f(X);
   fY = f(Y);
   fZ = f(Z);
   Image(:,:,1) = 116*fY - 16;    % L*
   Image(:,:,2) = 500*(fX - fY);  % a*
   Image(:,:,3) = 200*(fY - fZ);  % b*
end
return;

function Image = luv(Image,SrcSpace)
% Convert to CIE L*u*v* (CIELUV)
WhitePoint = [0.950456,1,1.088754];
WhitePointU = (4*WhitePoint(1))./(WhitePoint(1) + 15*WhitePoint(2) + 3*WhitePoint(3));
WhitePointV = (9*WhitePoint(2))./(WhitePoint(1) + 15*WhitePoint(2) + 3*WhitePoint(3));

Image = xyz(Image,SrcSpace); % Convert to XYZ
U = (4*Image(:,:,1))./(Image(:,:,1) + 15*Image(:,:,2) + 3*Image(:,:,3));
V = (9*Image(:,:,2))./(Image(:,:,1) + 15*Image(:,:,2) + 3*Image(:,:,3));
Y = Image(:,:,2)/WhitePoint(2);
L = 116*f(Y) - 16;
Image(:,:,1) = L;                        % L*
Image(:,:,2) = 13*L.*(U - WhitePointU);  % u*
Image(:,:,3) = 13*L.*(V - WhitePointV);  % v*
return;

function Image = lch(Image,SrcSpace)
% Convert to CIE L*ch
Image = lab(Image,SrcSpace);  % Convert to CIE L*ab
H = atan2(Image(:,:,3),Image(:,:,2));
H = H*180/pi + 360*(H < 0);
Image(:,:,2) = sqrt(Image(:,:,2).^2 + Image(:,:,3).^2);  % C
Image(:,:,3) = H;                                        % H
return;

function Image = huetorgb(m0,m2,H)
% Convert HSV or HSL hue to RGB
N = size(H);
H = min(max(H(:),0),360)/60;
m0 = m0(:);
m2 = m2(:);
F = H - round(H/2)*2;
M = [m0, m0 + (m2-m0).*abs(F), m2];
Num = length(m0);
j = [2 1 0;1 2 0;0 2 1;0 1 2;1 0 2;2 0 1;2 1 0]*Num;
k = floor(H) + 1;
Image = reshape([M(j(k,1)+(1:Num).'),M(j(k,2)+(1:Num).'),M(j(k,3)+(1:Num).')],[N,3]);
return;

function H = rgbtohue(Image)
% Convert RGB to HSV or HSL hue
[M,i] = sort(Image,3);
i = i(:,:,3);
Delta = M(:,:,3) - M(:,:,1);
Delta = Delta + (Delta == 0);
R = Image(:,:,1);
G = Image(:,:,2);
B = Image(:,:,3);
H = zeros(size(R));
k = (i == 1);
H(k) = (G(k) - B(k))./Delta(k);
k = (i == 2);
H(k) = 2 + (B(k) - R(k))./Delta(k);
k = (i == 3);
H(k) = 4 + (R(k) - G(k))./Delta(k);
H = 60*H + 360*(H < 0);
H(Delta == 0) = nan;
return;

function Rp = gammacorrection(R)
Rp = real(1.099*R.^0.45 - 0.099);
i = (R < 0.018);
Rp(i) = 4.5138*R(i);
return;

function R = invgammacorrection(Rp)
R = real(((Rp + 0.099)/1.099).^(1/0.45));
i = (R < 0.018);
R(i) = Rp(i)/4.5138;
return;

function fY = f(Y)
fY = real(Y.^(1/3));
i = (Y < 0.008856);
fY(i) = Y(i)*(841/108) + (4/29);
return;

function Y = invf(fY)
Y = fY.^3;
i = (Y < 0.008856);
Y(i) = (fY(i) - 4/29)*(108/841);
return;

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function h=vl_plotsiftdescriptor(d,f,varargin)
% VL_PLOTSIFTDESCRIPTOR   Plot SIFT descriptor
%   VL_PLOTSIFTDESCRIPTOR(D) plots the SIFT descriptor D. If D is a
%   matrix, it plots one descriptor per column. D has the same format
%   used by VL_SIFT().
%
%   VL_PLOTSIFTDESCRIPTOR(D,F) plots the SIFT descriptors warped to
%   the SIFT frames F, specified as columns of the matrix F. F has the
%   same format used by VL_SIFT().
%
%   H=VL_PLOTSIFTDESCRIPTOR(...) returns the handle H to the line
%   drawing representing the descriptors.
%
%   The function assumes that the SIFT descriptors use the standard
%   configuration of 4x4 spatial bins and 8 orientations bins. The
%   following parameters can be used to change this:
%
%   NumSpatialBins:: 4
%     Number of spatial bins in both spatial directions X and Y.
%
%   NumOrientationBins:: 8
%     Number of orientation bis.
%
%   MagnificationFactor:: 3
%     Magnification factor. The width of one bin is equal to the scale
%     of the keypoint F multiplied by this factor.
%
%   See also: VL_SIFT(), VL_PLOTFRAME(), VL_HELP().

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

opts.magnificationFactor = 3.0 ;
opts.numSpatialBins = 4 ;
opts.numOrientationBins = 8 ;
opts.maxValue = 0 ;

if nargin > 1
  if ~ isnumeric(f)
    error('F must be a numeric type (use [] to leave it unspecified)') ;
  end
end

opts = vl_argparse(opts, varargin) ;

% --------------------------------------------------------------------
%                                                  Check the arguments
% --------------------------------------------------------------------

if(size(d,1) ~= opts.numSpatialBins^2 * opts.numOrientationBins)
  error('The number of rows of D does not match the geometry of the descriptor') ;
end

if nargin > 1
  if (~isempty(f) & (size(f,1) < 2 | size(f,1) > 6))
    error('F must be either empty of have from 2 to six rows.');
  end

  if size(f,1) == 2
    % translation only
    f(3:6,:) = deal([10 0 0 10]') ;
    %f = [f; 10 * ones(1, size(f,2)) ; 0 * zeros(1, size(f,2))] ;
  end

  if size(f,1) == 3
    % translation and scale
    f(3:6,:) = [1 0 0 1]' * f(3,:) ;
    %f = [f; 0 * zeros(1, size(f,2))] ;
  end

  if size(f,1) == 4
    c = cos(f(4,:)) ;
    s = sin(f(4,:)) ;
    f(3:6,:) = bsxfun(@times, f(3,:), [c ; s ; -s ; c]) ;
  end

  if size(f,1) == 5
    assert(false) ;
    c = cos(f(4,:)) ;
    s = sin(f(4,:)) ;
    f(3:6,:) = bsxfun(@times, f(3,:), [c ; s ; -s ; c]) ;
  end

  if(~isempty(f) & size(f,2) ~= size(d,2))
    error('D and F have incompatible dimension') ;
  end
end

% Descriptors are often non-double numeric arrays
d = double(d) ;
K = size(d,2) ;

if nargin < 2 | isempty(f)
  f = repmat([0;0;1;0;0;1],1,K) ;
end

% --------------------------------------------------------------------
%                                                           Do the job
% --------------------------------------------------------------------

xall=[] ;
yall=[] ;

for k=1:K
  [x,y] = render_descr(d(:,k), opts.numSpatialBins, opts.numOrientationBins, opts.maxValue) ;
  xall = [xall opts.magnificationFactor*f(3,k)*x + opts.magnificationFactor*f(5,k)*y + f(1,k)] ;
  yall = [yall opts.magnificationFactor*f(4,k)*x + opts.magnificationFactor*f(6,k)*y + f(2,k)] ;
end

h=line(xall,yall) ;

% --------------------------------------------------------------------
function [x,y] = render_descr(d, numSpatialBins, numOrientationBins, maxValue)
% --------------------------------------------------------------------

% Get the coordinates of the lines of the SIFT grid; each bin has side 1
[x,y] = meshgrid(-numSpatialBins/2:numSpatialBins/2,-numSpatialBins/2:numSpatialBins/2) ;

% Get the corresponding bin centers
xc = x(1:end-1,1:end-1) + 0.5 ;
yc = y(1:end-1,1:end-1) + 0.5 ;

% Rescale the descriptor range so that the biggest peak fits inside the bin diagram
if maxValue
    d = 0.4 * d / maxValue ;
else
    d = 0.4 * d / max(d(:)+eps) ;
end

% We scramble the the centers to have them in row major order
% (descriptor convention).
xc = xc' ;
yc = yc' ;

% Each spatial bin contains a star with numOrientationBins tips
xc = repmat(xc(:)',numOrientationBins,1) ;
yc = repmat(yc(:)',numOrientationBins,1) ;

% Do the stars
th=linspace(0,2*pi,numOrientationBins+1) ;
th=th(1:end-1) ;
xd = repmat(cos(th), 1, numSpatialBins*numSpatialBins) ;
yd = repmat(sin(th), 1, numSpatialBins*numSpatialBins) ;
xd = xd .* d(:)' ;
yd = yd .* d(:)' ;

% Re-arrange in sequential order the lines to draw
nans = NaN * ones(1,numSpatialBins^2*numOrientationBins) ;
x1 = xc(:)' ;
y1 = yc(:)' ;
x2 = x1 + xd ;
y2 = y1 + yd ;
xstars = [x1;x2;nans] ;
ystars = [y1;y2;nans] ;

% Horizontal lines of the grid
nans = NaN * ones(1,numSpatialBins+1);
xh = [x(:,1)' ; x(:,end)' ; nans] ;
yh = [y(:,1)' ; y(:,end)' ; nans] ;

% Verical lines of the grid
xv = [x(1,:) ; x(end,:) ; nans] ;
yv = [y(1,:) ; y(end,:) ; nans] ;

x=[xstars(:)' xh(:)' xv(:)'] ;
y=[ystars(:)' yh(:)' yv(:)'] ;

author: yihui-he

repo: panorama

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


text:
function pdf = triangle_pdf ( x, a, b, c )

%*****************************************************************************80
%
%% TRIANGLE_PDF evaluates the Triangle PDF.
%
%  Discussion:
%
%    Given points A <= B <= C, the probability is 0 to the left of A,
%    rises linearly to a maximum of 2/(C-A) at B, drops linearly to zero
%    at C, and is zero for all values greater than C.
%
%  Formula:
%
%    PDF(A,B,C;X)
%      = 2 * ( X - A ) / ( B - A ) / ( C - A ) for A <= X <= B
%      = 2 * ( C - X ) / ( C - B ) / ( C - A ) for B <= X <= C.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%
%    Input, real A, B, C, the parameters of the PDF.
%    A <= B <= C and A < C.
%
%    Output, real PDF, the value of the PDF.
%
  if ( x <= a )

    pdf = 0.0;

  elseif ( x <= b )

    if ( a == b )
      pdf = 0.0;
    else
      pdf = 2.0 * ( x - a ) / ( b - a ) / ( c - a );
    end

  elseif ( x <= c )

    if ( b == c )
      pdf = 0.0;
    else
      pdf = 2.0 * ( c - x ) / ( c - b ) / ( c - a );
    end

  else
    pdf = 0.0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
filename='CantileverBeam_Triangle_Linear_Fine';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1 0.1];
constraint = {'volumeConstraint'};
optimizer = 'PROJECTED GRADIENT';
incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.5;
Perimeter_target=1;
optimality_final =1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function ind = sub2indn(varargin)

% function ind = sub2indn(siz,A,nanOutBound)
% ------------------------------------------------------------------------
% This function is similar to MATLAB's sub2ind function. However the input
% subscrip indices are provided as an nxm array, such that numel(siz)=m and
% numel(ind)=n. The output are the equivalent linear indices
%
%   See also: ind2subn ind2sub sub2ind.
%
% Change log
% 2018/05/17 Optionally output NaN values for out of bound indices
% ------------------------------------------------------------------------

%% Parse input

switch nargin
    case 2
        siz=varargin{1};
        A=varargin{2};
        nanOutBound=0;
    case 3
        siz=varargin{1};
        A=varargin{2};
        nanOutBound=varargin{3};
end

%%
siz = double(siz);
numDim = numel(siz);
k = cumprod(siz);

if numDim < 2
    error('Invalid size specified. The number of dimensions should be equal or larger than 2');
end

if any(size(A)==0)
    A=[];
end

if ~isempty(A)

    if size(A,2) ~= numDim
        error('The specified array size and number of subscript index columns do not match');
    end

    if ~nanOutBound
        if any(A(:)<1) || any(max(A,[],1)>siz)
            %Verify subscripts are within range
            error('Index out of range');
        end
    else
        A(A<1)=NaN;
        A(A>siz(ones(size(A,1),1),:))=NaN;
    end

    ind=A(:,1);
    for q=2:1:numDim
        ind = ind + (double(A(:,q))-1)*k(q-1);
    end
else
    ind=[];
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2019  Kevin Mattheus Moerman
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function theta=orderpoints(A, center, axis)
X0_in=center(1);
Y0_in=center(2);
x_axis=axis(:,1);
y_axis=axis(:,2);
centerlized=A-repmat([X0_in, Y0_in], size(A,1),1);
% length=sqrt(sum(centerlized.^2,2));
% cosx=centerlized*x_axis./length;
% cosy=centerlized*y_axis./length;
%
% for i=1:size(A,1)
%     if (cosy(i)>0)
%         theta(i)=acos(cosx(i));
%     else
%         theta(i)=2*pi-acos(cosx(i));
%
%     end
% end

x1=centerlized(:,1);
y1=centerlized(:,2);
x2=repmat(x_axis(1), numel(x1),1);
y2=repmat(x_axis(2), numel(x1),1);
theta = mod(atan2(x1.*y2-x2.*y1,x1.*x2+y1.*y2),2*pi)';

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [gout,a,fc,L,info] = waveletfilters(Ls, scales, varargin)
%WAVELETFILTERS Generates wavelet filters
%   Usage: [gout,a,fc,L,info] = waveletfilters(Ls,scales)
%          [gout,a,fc,L,info] = waveletfilters(Ls,'bins', fs, fmin, fmax, bins)
%          [gout,a,fc,L,info] = waveletfilters(Ls,'linear', fs, fmin, fmax, channels)
%
%   Input parameters:
%         Ls     : System length
%         scales : Vector of wavelet scales
%   Output parameters:
%         gout  : Cell arrary of wavelet filters
%         a     : Downsampling rate for each channel.
%         fc    : Center frequency of each channel.
%         L     : Next admissible length suitable for the generated filters.
%         info  : Struct with additional outputs
%
%   `waveletfilters(Ls,scales)` constructs a system of wavelet filters covering
%   scales in the range *scales* for system length *Ls*. A scale of 1 corresponds
%   to a wavelet filter with peak positioned at the frequency 0.1 relative
%   to the Nyquist rate.
%
%   `[g,a,fc]=waveletfilters(Ls, 'bins', fs,fmin,fmax,bins)` constructs a set of
%   wavelets *g* which cover the required frequency range
%   `fmin`-`fmax` with `bins` filters per octave starting at `fmin`. All
%   filters have (approximately) equal $Q=f_c/f_b$. The
%   frequency interval below fmin not covered by these is captured
%   by additional lowpass filter(s) . The signal length *Ls*
%   is mandatory, since we need to avoid too narrow frequency windows
%
%   `[g,a,fc]=waveletfilters(Ls, 'linear', fs, fmin, fmax, channels)` constructs
%    a set of wavelets *g* which cover the required frequency range
%   `fmin`-`fmax` with `channels` equidistantly spaced filters starting at `fmin`.
%
%   Wavelet types
%   --------------------
%
%   The following wavelets can be passed as a flag:
%
%   'cauchy'     Cauchy wavelet (default parameters [alpha beta gamma] = [300 0 3])
%
%   'morse'      Generalized morse wavelet (default parameters [alpha beta gamma] = [300 0 3])
%
%   'morlet'     Morlet wavelet (default parameters sigma = [4])
%
%   'fbsp'       Frequency B-spline wavelet (default parameters [order fb] = [4 2])
%
%   'analyticsp' Analytic spline wavelet (default parameters [order fb] = [4 2])
%
%   'cplxsp'     Complex spline wavelet (default parameters [order fb] = [4 2])
%
%   A scale of 1 corresponds
%   to a wavelet filter with peak positioned at the frequency 0.1 relative
%   to the Nyquist rate. By default, this function does not allow center frequencies
%   exceeding the Nyquist rate, except with the optional parameter 'analytic',
%   see below. This implies that scales below 0.1 (or 0.05 for the 'analytic' scheme)
%   are not supported.
%   For more details on the construction of the wavelets and the available
%   wavelet types, please see |freqwavelet|.
%
%   By default, wavelet filters are peak normalized before being adjusted
%   to the proposed downsampling factor. The peak normalization can be
%   overridden by forwarding any norm flag accepted by |setnorm|.
%
%   Downsampling factors
%   --------------------
%
%   The integer downsampling rates of the channels must all divide the
%   signal length, |filterbank| will only work for input signal lengths
%   being multiples of the least common multiple of the downsampling rates.
%   See the help of |filterbanklength|.
%   The fractional downsampling rates restrict the filterbank to a single
%   length *L=Ls*.
%
%   `[gout,a]=waveletfilters(...,'regsampling')` constructs a non-uniform
%   filterbank with integer subsampling factors. This is the default.
%
%   `[gout,a]=waveletfilters(...,'uniform')` constructs a uniform filterbank
%   where the integer downsampling rate is the same for all the channels. This
%   results in the most redundant representation which produces nice plots.
%
%   `[gout,a]=waveletfilters(...,'fractional')` constructs a filterbank with
%   fractional downsampling rates *a*.
%   This results in the least redundant system.
%
%   `[gout,a]=waveletfilters(...,'fractionaluniform')` constructs a filterbank
%   with fractional downsampling rates *a*, which are uniform for all filters
%   except the "filling" low-pass filter which can have different
%   fractional downsampling rates. This is useful when uniform subsampling
%   and low redundancy at the same time are desirable.
%
%   Lowpass filters
%   --------------------
%
%   `[gout,a]=waveletfilters(...,'single')` uses a single lowpass filter
%   for covering the range from zero frequency to the center frequency of the
%   largest scale specified. This is the default.
%
%   `[gout,a]=waveletfilters(...,'repeat')` constructs frequency-shifted
%   copies of the largest scale wavelet to cover the range from zero frequency
%   to the center frequency of the largest scale specified.
%
%   `[gout,a]=waveletfilters(...,'none')` foregoes the construction of a
%   lowpass filter. This option cannot be expected to yield an invertible
%   filterbank.
%
%   Additional parameters
%   ---------------------
%
%   `waveletfilter` accepts the following optional parameters:
%
%     'redmul',redmul    Redundancy multiplier. Increasing the value of this
%                        will make the system more redundant by lowering the
%                        channel downsampling rates. It is only used if the
%                        filterbank is a non-uniform filterbank. Default
%                        value is *1*. If the value is less than one, the
%                        system may no longer be painless.
%
%     'redtar',redtar    Target redundancy. The downsampling factors will be
%                        adjusted to achieve a redundancy as close as possible
%                        to 'redtar'.
%
%     'trunc_at',trunc_at     Applies hard thresholding of the wavelet filters
%                             at the specified threshold value to reduce their
%                             support size.
%                             The default value is *trunc_at=10e-5*. When no
%                             truncation is desired, *trunc_at=0* should be chosen.
%
%     'delay',delay      A scalar, numeric vector of function handle that
%                        specifies delays for the wavelet filters. A
%                        numeric vector must have at least as many entries
%                        as there are filters in the filterbank. A function
%                        handle must accept two inputs *(k-1,a(k))*, where
%                        *k* is the channel index and *a* are the
%                        downsampling rates. If a function handle is given
%                        and 'redtar' is specified, delays are computed
%                        based on the final value of *a*.
%
%     'real'             Allows positive scales with center frequencies up
%                        to Nyquist. This is the default.
%
%     'complex'          Allows positive scales with center frequencies up
%                        to Nyquist, which are also mirrored to cover
%                        negative scales.
%
%     'analytic'         Allows positive scales with center frequencies up
%                        to twice the Nyquist frequency. This setting is
%                        suitable for the analysis of analytic signals.
%
%     'startfreq'        Allows to manually set a starting frequency for
%                        the wavelet range. Can not be lower than fmin.
%
%   Examples:
%   ---------
%
%   In the first example, we analyze a glockenspiel signal with a
%   regularly sampled wavelet filterbank using a frequency B-spline
%   wavelet of order 4 and with parameter fb=3 and visualize the result:::
%
%     [f,fs]=gspi;  % Get the test signal
%     Ls = length(f);
%     scales = linspace(10,0.1,100);
%     [g,a,fc,L]=waveletfilters(Ls,scales, {'fbsp', 4, 3}, 'repeat');
%     c=filterbank(f,g,a);
%     plotfilterbank(c,a,fc,fs,90);
%
%   In the second example, we construct a wavelet filterbank with a
%   lowpass channels based on a Cauchy wavelet and verify it.
%   The plot shows the frequency responses of
%   filters used for analysis (top) and synthesis (bottom). :::
%
%     [f,fs]=greasy;  % Get the test signal
%     Ls = length(f);
%     M0 = 511; %Desired number of channels (without 0 Hz-lowpass channel)
%     max_freqDiv10 = 10;  % 10 corresponds to the nyquist frequency
%     freq_step = max_freqDiv10/M0;
%     rate = 44100;
%     start_index = 1;
%     min_freqHz = rate/10*freq_step
%     min_scale_freq = min_freqHz*start_index
%     min_freqDiv10 = freq_step*start_index; %1/25; % By default, the reference scale for freqwavelet has center frequency 0.1
%     scales = 1./linspace(min_freqDiv10,max_freqDiv10,M0);
%     alpha = 1-2/(1+sqrt(5)); % 1-1/(goldenratio) delay sequence
%     delays = @(n,a) a*(mod(n*alpha+.5,1)-.5);
%     CauchyAlpha = 600;
%     [g, a,fc,L,info] = waveletfilters(Ls,scales,{'cauchy',CauchyAlpha},'uniform','single','energy', 'delay',delays, 'redtar', 8);
%
%     c=filterbank(f,{'realdual',g},a);
%     r=2*real(ifilterbank(c,g,a));
%     if length(r) > length(f)
%         norm(r(1:length(f))-f)
%     else
%         norm(r-f(1:length(r)))
%      end
%     % Plot frequency responses of individual filters
%     gd=filterbankrealdual(g,a,L);
%     figure(1);
%     subplot(2,1,1);
%     filterbankfreqz(gd,a,L,fs,'plot','linabs','posfreq');
%
%     subplot(2,1,2);
%     filterbankfreqz(g,a,L,fs,'plot','linabs','posfreq');
%
%   See also: freqwavelet, filterbank, setnorm

% AUTHORS: Nicki Holighaus, Zdenek Prusa, Guenther Koliander, Clara Hollomey

complainif_notenoughargs(nargin,2,upper(mfilename));
complainif_notposint(Ls,'Ls',upper(mfilename));

%parse input arguments:
if ~isnumeric(scales)
    fs = varargin{1};
    fmin = varargin{2};
    fmax = varargin{3};
    channels = varargin{4};
    switch scales
        case 'linear'
            definput.flags.inputmode = {'linear', 'logarithmic', 'bins', 'scales'};
        %case 'logarithmic'
        %    definput.flags.inputmode = {'logarithmic', 'bins', 'scales', 'linear'};
        case 'bins'
            definput.flags.inputmode = {'bins', 'scales', 'linear', 'logarithmic'};
        otherwise
            error('%s: second argument must either be a scales vector or define the f-mapping.',upper(mfilename))
    end
    %this is a slightly more efficient way to remove the first 4 args from
    %varargin
    varargin = circshift(varargin,-4);
    varargin = varargin(1:end-4);
else
    definput.flags.inputmode = {'scales', 'linear', 'logarithmic', 'bins'};
end

definput.import={'setnorm'};
definput.importdefaults={'null'};
definput.flags.real = {'real','complex','analytic'};
definput.flags.lowpass  = {'single','repeat','none'};
definput.flags.sampling = {'regsampling','uniform',...
                           'fractional','fractionaluniform'};
definput.flags.wavelettype = getfield(arg_freqwavelet(),'flags','wavelettype');
definput.keyvals.redmul=1;
definput.keyvals.redtar=[];
definput.keyvals.delay = 0;
definput.keyvals.trunc_at  = 10^(-5);
definput.keyvals.fs = 2;
definput.keyvals.startfreq = [];%only relevant if 'repeat'

[varargin,winCell] = arghelper_filterswinparser(definput.flags.wavelettype,varargin);
[flags,kv]=ltfatarghelper({},definput,varargin);

if isempty(winCell), winCell = {flags.wavelettype}; end

if ~isa(kv.delay,'function_handle') && ~isnumeric(kv.delay)
    error('%s: delay must be a function handle or numeric.',upper(mfilename));
end

if ~isscalar(kv.redmul) || kv.redmul <= 0
    error('%s: redmul must be a positive scalar.',upper(mfilename));
end

if ~isempty(kv.redtar)
    if ~isscalar(kv.redtar) || kv.redtar <= 0
        error('%s: redtar must be a positive scalar.',upper(mfilename));
    end
end

%parse the input format: map fmin and fmax to scales according to the input
%parameter specification
if ~flags.do_scales
    nf = fs/2;
    if flags.do_linear
        min_freq = fmin/nf *10;%map to freqwavelets nyquist f
        max_freq = fmax/nf * 10;
        scales = 1./linspace(min_freq,max_freq,channels);
  % elseif flags.do_logarithmic
%
%        fc = 2.^linspace(log2(fmin), log2(fmax), channels);
%        fc = fc/nf * 10;
%        scales = 1./fc;

    elseif flags.do_bins

        if isscalar(channels)
            % Number of octaves
            b = ceil(log2(fmax/fmin))+1;
            bins = channels*ones(b,1);
        else
            bins = channels;
        end

        fc = zeros(sum(bins),1);

        ll = 0;
        for kk = 1:length(bins)
            fc(ll+(1:bins(kk))) = ...
                fmin*2.^(((kk-1)*bins(kk):(kk*bins(kk)-1)).'/bins(kk));
            ll = ll+bins(kk);
        end

        % Get rid of filters with frequency centers >=fmax and nf
        % This will leave the first bigger than fmax it it is lower than nf
        temp = find(fc>=fmax ,1);
        if fc(temp) >= nf
            fc = fc(1:temp-1);
        else
            fc = fc(1:temp);
        end

        channels = length(fc);
        min_freq = fmin/nf *10;%map to freqwavelets nyquist f
        max_freq = fmax/nf * 10;
        scales = 1./linspace(min_freq,max_freq,channels);
    end
    scales_sorted = sort(scales,'descend');
    if ~isempty(kv.startfreq)%set the start frequency
        startfreq = kv.startfreq/nf * 10;
        scales_start = find(1./scales_sorted > startfreq,1,'first');%find first scale whose equiv. f is larger than fmin
        scales = scales(scales_start:end);
    end
end

if ~isnumeric(scales) || any(scales < 0.1)
   error('%s: scales must be positive and numeric.',upper(mfilename));
end

if size(scales,2)>1
    if size(scales,1)==1
        % scales was a row vector.
        scales=scales(:);
    else
        error('%s: scales must be a vector.',upper(mfilename));
    end
end

%% Generate mother wavelet to determine parameters from
[~,info] = freqwavelet(winCell,Ls,1,'asfreqfilter','efsuppthr',kv.trunc_at,'basefc',0.1);
basea = info.aprecise;

%% Determine total number of filters and natural subsampling factor for lowpass
%[aprecise, M, lowpass_number, lowpass_at_zero] = c_det_lowpass(Ls, scales, basea, flags, kv);
if numel(scales) < 4 && flags.do_single
    error('%s: Lowpass generation requires at least 4 scales.',upper(mfilename));
elseif numel(scales) < 2 && flags.do_repeat
    error('%s: Lowpass generation requires at least two scales.',upper(mfilename));
end

% Get number of scales and sort them
M = numel(scales);
scales_sorted = sort(scales,'descend');
%% Determine total number of filters and natural subsampling factor for lowpass
if flags.do_repeat
% Maybe adjust this to not guarantee some distance between first filter and zero frequency.
    lowpass_number = scales_sorted(2)/(scales_sorted(1)-scales_sorted(2));
        if abs(lowpass_number - round(lowpass_number)) < eps*10^3
            % determine if lowpass is centered around 0 Hz
            lowpass_number = round(lowpass_number);
            lowpass_at_zero = 1;
        else
            lowpass_at_zero = 0;
        end
        lowpass_number = floor(lowpass_number);
        if lowpass_number == 0
            lowpass_number = 1;
        end
        M = M + lowpass_number;
        aprecise = (basea.*scales_sorted(1))*ones(lowpass_number,1);

elseif flags.do_single
    lowpass_number = 1;
    lowpass_at_zero = 1;
    M = M+1;
    %this is an estimated value
    aprecise = (0.2./scales_sorted(4))*Ls; % This depends on how single lowpass is called (l.195ff). Maybe automate. Adapt if necessary!!!
else
    lowpass_number = 0;
    lowpass_at_zero = 0;
    aprecise = [];
end

%% Get subsampling factors
aprecise = [aprecise;basea.*scales];

if any(aprecise<1)
    error(['%s: Bandwidth of at least one of the filters is bigger than fs. '],upper(mfilename));
end

aprecise=aprecise/kv.redmul;
if any(aprecise<1)
    error('%s: The maximum redundancy mult. for this setting is %5.2f',...
         upper(mfilename), min(basea./scales));
end

%% Compute the downsampling rate
if flags.do_regsampling
    a = ones(M,1);

    [lower_scale,~] = max(scales);
    [upper_scale,~] = min(scales);
    lower_scale = floor(log2(1/lower_scale));
    upper_scale = floor(log2(1/upper_scale));

    % Find minimum a in each octave and floor23 it
    % to shrink "a" to the next composite number
    ct=1;
    for kk = lower_scale:upper_scale
        tempidx = find( floor(log2(1./scales)) == kk );
        [~,tempminidx] = min(1/scales(tempidx));
        idx = tempidx(tempminidx);

        % Deal the integer subsampling factors
        a(tempidx) = floor23(aprecise(idx));
        ct=ct+1;
    end

    % Determine the minimal transform length lcm(a)
    L = filterbanklength(Ls,a);

    % Heuristic trying to reduce lcm(a)
    while L>2*Ls && ~(all(a==a(1)))
        maxa = max(a);
        a(a==maxa) = 0;
        a(a==0) = max(a);
        L = filterbanklength(Ls,a);
    end

elseif flags.do_fractional
    L = Ls;
    N=ceil(Ls./aprecise);
    a=[repmat(Ls,M,1),N];
elseif flags.do_fractionaluniform
    L = Ls;
    if lowpass_at_zero
        aprecise(2:end)= min(aprecise(2:end));
    else
        aprecise= repmat(min(aprecise),numel(aprecise),1);
    end
    N=ceil(Ls./aprecise);
    a=[repmat(Ls,M,1),N];
elseif flags.do_uniform
    a=floor(min(aprecise));
    L=filterbanklength(Ls,a);
    a = repmat(a,M,1);
end
% Get an expanded "a" / Convert "a" to LTFAT 2-column fractional format
afull=comp_filterbank_a(a,M,struct());

%if flags.do_uniform
%    a = a(:,1);
%end
%==========================================================================
%% Adjust the downsampling rates in order to achieve 'redtar'

if ~isempty(kv.redtar)
   if size(afull,2) == 2
        a = afull(:,1)./afull(:,2);
   else
       a = afull;
   end

    if ~flags.do_real
        org_red = sum(1./a);
    elseif lowpass_at_zero
        org_red = 1./a(1) + sum(2./a(2:end));
    else
        org_red = sum(2./a);
    end

    a = floor(a*org_red/kv.redtar);
    a(a==0) = 1;

    if ~flags.do_uniform
        N_new=ceil(L./a);
        if flags.do_complex
            N_new = [N_new;N_new(end:-1:2)];
        end
        a=[repmat(L,numel(N_new),1),N_new];
    else
        L = filterbanklength(L,a);
        a=[a,ones(length(a), 1)];
    end
else
    a = afull;
end

%% Compute the scaling of the filters and the numeric delay vector
% Filters are scaled such that the energy of the subband coefficients
% remains approximately constant independent of the decimation factor
if isa(kv.delay,'function_handle')
    delayvec = zeros(M,1);
    for kk = 1:M
        delayvec(kk) = kv.delay(kk-1,a(kk,1)./a(kk,2));
    end
elseif numel(kv.delay) == 1
    delayvec = repmat(kv.delay,M,1);
elseif ~isempty(kv.delay) && size(kv.delay,2) > 1
    delayvec = kv.delay(:);
else
    error('%s: delay must be scaler or have enough elements to cover all channels.',upper(mfilename));
end
scal=sqrt(a(:,1)./a(:,2));

if flags.do_complex

    if lowpass_at_zero
       a=[a;flipud(a(2:end,:))];
       scal=[scal;flipud(scal(2:end))];
       delayvec=[delayvec;flipud(delayvec(2:end))];
    else
        a=[a;flipud(a)];
        scal=[scal;flipud(scal)];
        delayvec=[delayvec;flipud(delayvec)];
    end

    [gout_positive,info_positive] = freqwavelet(winCell,L,scales,...
        'asfreqfilter','efsuppthr',kv.trunc_at,'basefc',0.1,...
        'scal',scal(lowpass_number+1:M),'delay', delayvec(lowpass_number+1:M),flags.norm);
    [gout_negative,info_negative] = freqwavelet(winCell,L,-flipud(scales),...
        'asfreqfilter','efsuppthr',kv.trunc_at,'basefc',0.1,...
        'negative','scal',scal(M+1:M+numel(scales)),'delay', delayvec(M+1:M+numel(scales)), flags.norm);
    gout = [gout_positive,gout_negative];
    fields = fieldnames(info_positive);
    info = struct();
    for kk = 1:length(fields)
            info.(fields{kk}) = [info_positive.(fields{kk}),info_negative.(fields{kk})];
    end
elseif flags.do_analytic
    [gout,info] = freqwavelet(winCell,L,scales,...
        'asfreqfilter','efsuppthr',kv.trunc_at,'basefc',0.1,...
        'analytic','scal',scal(lowpass_number+1:M),'delay', delayvec(lowpass_number+1:M),flags.norm);
else
    if lowpass_at_zero
        % Scale the lowpass filters
        scal(1)=scal(1)/sqrt(2);
    end

    [gout,info] = freqwavelet(winCell,L,scales,'asfreqfilter','efsuppthr',...
        kv.trunc_at,'basefc',0.1,'scal',scal(lowpass_number+1:M),'delay', delayvec(lowpass_number+1:M),flags.norm);
end

%% Generate lowpass filters if desired
[gout, info] = comp_fblowpassfilters(winCell, gout, a, L, info, scales, scal, delayvec(1:lowpass_number), lowpass_at_zero, kv, flags);

info.lowpassstart = lowpass_number + 1;%startindex of actual wavelets (tentative)
% Assign fc and adjust for sampling rate
if flags.do_scales
    fc = (kv.fs/2).*info.fc;
else
    fc = nf.*info.fc;
end

if flags.do_uniform || flags.do_regsampling
    a = a(:,1);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [opts, args] = vl_argparse(opts, args, varargin)
%VL_ARGPARSE Parse list of parameter-value pairs.
%   OPTS = VL_ARGPARSE(OPTS, ARGS) updates the structure OPTS based on
%   the specified parameter-value pairs ARGS={PAR1, VAL1, ... PARN,
%   VALN}. If a parameter PAR cannot be matched to any of the fields
%   in OPTS, the function generates an error.
%
%   Parameters that have a struct value in OPTS are processed
%   recursively, updating the individual subfields.  This behaviour
%   can be suppressed by using VL_ARGPARSE(OPTS, ARGS, 'nonrecursive'),
%   in which case the struct value is copied directly (hence deleting any
%   existing subfield existing in OPTS). A direct copy occurrs also if the
%   struct value in OPTS is a structure with no fields. The nonrecursive
%   mode is especially useful when a processing time is a concern.
%
%   One or more of the (PAR, VAL) pairs in the argument list can be
%   replaced by a structure; in this case, the fields of the structure
%   are used as paramater names and the field values as parameter
%   values. This behaviour, while orthogonal to structure-valued parameters,
%   is also disabled in the 'nonrecursive' mode.
%
%   [OPTS, ARGS] = VL_ARGPARSE(OPTS, ARGS) copies any parameter in
%   ARGS that does not match OPTS back to ARGS instead of producing an
%   error. Options specified as structures are passed back as a list
%   of (PAR, VAL) pairs.
%
%   Example::
%     The function can be used to parse a list of arguments
%     passed to a MATLAB functions:
%
%        function myFunction(x,y,z,varargin)
%        opts.parameterName = defaultValue ;
%        opts = vl_argparse(opts, varargin)
%
%     If only a subset of the options should be parsed, for example
%     because the other options are interpreted by a subroutine, then
%     use the form
%
%        [opts, varargin] = vl_argparse(opts, varargin)
%
%     that copies back to VARARGIN any unknown parameter.
%
%   See also: VL_HELP().

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

if ~isstruct(opts) && ~isobject(opts), error('OPTS must be a structure') ; end
if ~iscell(args), args = {args} ; end

recursive = true ;
if numel(varargin) == 1
  if strcmp(lower(varargin{1}), 'nonrecursive') ;
    recursive = false ;
  else
    error('Unknown option specified.') ;
  end
end
if numel(varargin) > 1
  error('There can be at most one option.') ;
end

optNames = fieldnames(opts)' ;

% convert ARGS into a structure
ai = 1 ;
keep = false(size(args)) ;
while ai <= numel(args)

  % Check whether the argument is a (param,value) pair or a structure.
  if recursive && isstruct(args{ai})
    params = fieldnames(args{ai})' ;
    values = struct2cell(args{ai})' ;
    if nargout == 1
      opts = vl_argparse(opts, vertcat(params,values)) ;
    else
      [opts, rest] = vl_argparse(opts, reshape(vertcat(params,values), 1, [])) ;
      args{ai} = cell2struct(rest(2:2:end), rest(1:2:end), 2) ;
      keep(ai) = true ;
    end
    ai = ai + 1 ;
    continue ;
  end

  if ~isstr(args{ai})
    error('Expected either a param-value pair or a structure.') ;
  end

  param = args{ai} ;
  value = args{ai+1} ;

  p = find(strcmpi(param, optNames)) ;
  if numel(p) ~= 1
    if nargout == 1
      error('Unknown parameter ''%s''', param) ;
    else
      keep([ai,ai+1]) = true ;
      ai = ai + 2 ;
      continue ;
    end
  end
  field = optNames{p} ;

  if ~recursive
    opts.(field) = value ;
  else
    if isstruct(opts.(field)) && numel(fieldnames(opts.(field))) > 0
      % The parameter has a  non-empty struct value in OPTS:
      % process recursively.
      if ~isstruct(value)
        error('Cannot assign a non-struct value to the struct parameter ''%s''.', ...
          field) ;
      end
      if nargout > 1
        [opts.(field), args{ai+1}] = vl_argparse(opts.(field), value) ;
      else
        opts.(field) = vl_argparse(opts.(field), value) ;
      end
    else
      % The parameter does not have a struct value in OPTS: copy as is.
      opts.(field) = value ;
    end
  end

  ai = ai + 2 ;
end

args = args(keep) ;

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function [xr,xc,xm] = z2x(z,lc,lm)
if length(z)==12
    xr = [z(1:2),0,0,0,z(3)]; a = z(4:7); th = z(8:11); w = z(12);
else
    xr = [z(1:2),0,0,0,0]; a = z(3)*ones(1,4); th = zeros(1,4); w = z(4);
end
xc = xr2c(xr,w);
xm = xr2m(xr,a,th,lc,lm);
end

repo: paper-simulation


text:
function [] = Resultant_Oscillation( Nb,n )

% Nb=10;
% n=10;
% figure(1)
% hold on

R=100;

psi=[0:1:360]*pi/180;

npsi=length(psi);

bih1=sin(psi);

for i=0:Nb-1
    bihn(i+1,:)=sin(n*(2*i*pi/Nb+psi));
end
resultant=sum(bihn);
if max(resultant)>1e-3
    resultant=resultant/max(resultant);
end

% title('Resultant Oscillation Transmitted to Hull','Position',[179.129 1.007 17.321],'FontWeight','Bold')
hold on
axis([0,360,-1,1])
Waves=plot(psi*180/pi,resultant,'r',  psi*180/pi,bih1,'b');
% legend('Transmitted','1 Per rev','Location','BestOutside')
set(gca,'ytick',0)
set(gca,'xtick',[0:90:360])
set(gca,'XAxisLocation','top')
set(gca,'XTickLabel',[])
set(gca,'YTickLabel',[])

grid on
save DataFile1 Waves


text:
function vw = computeContrastMap2(vw, active, control, contrastName, varargin);
%
%  view = computeContrastMap2(view, active, control, contrastName, [options]);
%
% Compute a statistical contrast map on data for a view, applying
% a GLM if necessary. Updated version, using the rewritten GLM code.
%
% The map is computed off data from the scan group assigned to
% the view's current scan. (You can specify a different scan
% group as an optional argument). The map is added to the
% view's map field, and saved in the view's data dir under
% the specified contrast name.
%
% By default, the map represents the -log10(p) value of a statistical
% T test between the active and control conditions. So, a value of
% 2 indicates that the active conditions are greater than the
% control with a significance of 10^-2, a value of 3 is 10^-3, and
% a value of -2 indicates that control > active with a likelihood of
% 10^-2. Other tests (F test rather than T test) and mapUnits are possible
% using the options specified below.
%
%
%
%
% Options include:
%
%   [add options for varying slice, statistical test, mapUnits]
%
%
% ras 06/06
% ras 02/07: had some 'quickCompute' code in here, but after testing out
% load times, it looks like loading part of the GLM information doesn't
% save time over loading everything. So, kept it simple and removed the
% quickCompute code. Also: saves the mapUnits in a separate variable, will
% update param map code to show mapUnits if they're specified (and maybe if
% the user sets a flag).
% dar 02/07: made the bicolormap assignment conditional on the presence of
% a ui field in view (was crashing on hiddeninplanes).
if nargin<1, help(mfilename); return; end

%%%%% default params
test = 't';

% allow for a user input if the active and control
% conditions aren't specified
if ~exist('active','var') | isempty(active) |  ...
	~exist('control','var') | isempty(control)
   prompt={'Active (+) Conditions:' ...
           'Control (-) Conditions:' ...
           'Contrast Name (if blank, will name based on condition names)'};
   def={'1' '0' ''};
   dlgTitle='Compute Contrast Map 2';

   answer=inputdlg(prompt, dlgTitle, 1, def);

   active = str2num(answer{1});
   control = str2num(answer{2});
   contrastName = answer{3};
end
scan = viewGet(vw,'curScan');
opts = {};
verbose = prefsVerboseCheck;
if verbose,		tic;	end   % time the map computation
tcWeights = [];

%%%%% parse the options
for i = 1:length(varargin)
    switch(lower(varargin{i}))
        case 'f',           test = 'f';
        case 't',           test = 't';
		case 'fdr',		opts = [opts {'mapUnits' 'fdr'}];
        case {'w' 'weights'}  % contrast weights
            opts = [ opts { 'weights' varargin{i+1} } ];
        case {'tcweights'}, tcWeights = varargin{i+1};
        case 'mapunits', opts = [opts {'mapUnits'}];, opts=[opts {varargin{i+1}}];
    end
end

%% auto-add tcWeights option for 'selective averaging' GLMs:
params = er_getParams(vw);
if params.glmHRF==0  % flag for deconvolution
	% what this means is, if the GLM didn't assume an HRF, but instead
	% generated several beta values for each condition (1 per time point in a
	% peri-stimulus time window), we want to specify those betas (=frames in time
	% course) to use for the contrast. This should be the same points selected
	% as the 'peakPeriod' parameter. Formally, tcWeights should support
	% differential weighting across time ponts, but right now the glm_contrast code
	% treats it as a list of frames:
	f1 = fix( min(params.timeWindow) / params.framePeriod );
	f2 = fix( max(params.timeWindow) / params.framePeriod );
	pk1 = fix( min(params.peakPeriod) / params.framePeriod );
	pk2 = fix( max(params.peakPeriod) / params.framePeriod );
	tcWeights = find(ismember(f1:f2, pk1:pk2));
	opts = [opts {'tcWeights' tcWeights}];
end

%% default contrast name
if ~exist('contrastName','var') | isempty(contrastName)
	trials = er_concatParfiles(vw);
    activeNames = ''; controlNames = '';
    for i = 1:length(active)
        j = find(trials.condNums==active(i));
        activeNames = [activeNames trials.condNames{j}];
    end
    for i = 1:length(control)
        j = find(trials.condNums==control(i));
        controlNames = [controlNames trials.condNames{j}];
    end
    contrastName = sprintf('%sV%s', activeNames, controlNames)
end

%% check if a GLM has been run; if not, run one:
checkPath = fullfile(dataDir(vw),sprintf('Scan%i',scan),'glmSlice1.mat');
if ~exist(checkPath,'file')
    % no GLM file found -- ask user if they want to run one
    mrGlobals;
    [scans dt] = er_getScanGroup(vw,scan);
    names = {dataTYPES.name};
    q = 'A GLM wasn''t found for this slice and scan. ';
    q = [q 'Do you want to run one now? '];
    q = [q sprintf('%s scans %s',names{dt},num2str(scans))];
    resp = questdlg(q,'Compute Contrast Map');
    if isequal(resp,'Yes')
        % compute one
        vw = selectDataType(vw,dt);
        vw = setCurScan(vw,scans(1));
        [newDt newScan] = applyGlm(vw, scans);
        vw  = selectDataType(vw, newDt);
        vw = setCurScan(vw, newScan);
    else
        disp('Aborting computeContrastMap2')
        return
    end
end

% initialize map cell (of size nScans):
% If another map with the same name has already been
% computed, load this: we'll just plug in the data for
% the current scan:
mapPath = fullfile(dataDir(vw),[contrastName '.mat']);
if exist(mapPath, 'file')
    load(mapPath, 'map', 'mapName');
else
    mapName = contrastName;
    map = cell(1, viewGet(vw, 'numScans'));
end

% initalize volume for the map
mapVol = NaN*ones(dataSize(vw, scan));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% main loop: load data from each slice and compute map  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if verbose, hwait = mrvWaitbar(0,'Computing Contrast Map'); end

nSlices = viewGet(vw, 'numSlices');
for slice = 1:nSlices
	model  = loadGlmSlice(vw, slice, scan);

	mapVol(:,:,slice) = glm_contrast(model, active, control, ...
									   'size', viewGet(vw, 'sliceDims'), opts);

	if verbose, mrvWaitbar(slice/nSlices, hwait);  end
end

if verbose, close(hwait);	end

%% set new map in view
mapUnits = '-log(p)';  % default
if cellfind(opts, 'mapUnits')
    mapUnits = opts{ cellfind(opts, 'mapUnits') + 1 };
end
if isequal(mapUnits, 'log10p'), mapUnits = '-log(p)'; end   % clarify
map{scan} = mapVol;
vw.mapUnits = mapUnits;
vw = setParameterMap(vw, map, mapName);
if checkfields(vw, 'ui', 'mapMode')
	vw = bicolorCmap(vw);  % take a chance that people will like this
end

% set the variance explained as the coherence
varExpPath = fullfile(dataDir(vw), 'Proportion Variance Explained');
if check4File(varExpPath)
	tmp = load(varExpPath, 'map');
	co = tmp.map;
else
	% this will wipe out any existing co data when the map is loaded:
	% but since we're in the GLMs data type, that should be ok.
	co = cell(1, numScans(vw));
end

%% save results
if exist(mapPath,'file')
    save(mapPath, 'map', 'mapName', 'mapUnits', 'co', '-append');
else
    save(mapPath, 'map', 'mapName', 'co', 'mapUnits');
end
fprintf('Saved Contrast Map in %s.\n', mapPath);

%% refresh screen according to whether we're using mrVista 1 or 2:
global GUI
if ~isempty(GUI)
    sessionGUI_selectDataType;  % mrVista 2
else
    vw = refreshScreen(vw, 1);        % mrVista 1
end

if verbose
	fprintf('Time to compute map: %i min, %3.1f sec.\n\n', ...
		floor(toc/60), mod(toc, 60));
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% eegplugin_firfilt() - EEGLAB plugin for filtering data using linear-
%                       phase FIR filters
%
% Usage:
%   >> eegplugin_firfilt(fig, trystrs, catchstrs);
%
% Inputs:
%   fig        - [integer]  EEGLAB figure
%   trystrs    - [struct] "try" strings for menu callbacks.
%   catchstrs  - [struct] "catch" strings for menu callbacks.
%
% Author: Andreas Widmann, University of Leipzig, Germany, 2005

%123456789012345678901234567890123456789012345678901234567890123456789012

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

function vers = eegplugin_firfilt(fig, trystrs, catchstrs)

    vers = 'firfilt1.6.1';
    if nargin < 3
        error('eegplugin_firfilt requires 3 arguments');
    end

    % add folder to path
    % -----------------------
    if ~exist('pop_firws')
        p = which('eegplugin_firfilt');
        p = p(1:findstr(p,'eegplugin_firfilt.m')-1);
        addpath([p vers]);
    end

    % find import data menu
    % ---------------------
    menu = findobj(fig, 'tag', 'filter');

    % menu callbacks
    % --------------
    comfirfiltnew = [trystrs.no_check '[EEG LASTCOM] = pop_eegfiltnew(EEG);' catchstrs.new_and_hist];
    comfirws = [trystrs.no_check '[EEG LASTCOM] = pop_firws(EEG);' catchstrs.new_and_hist];
    comfirpm = [trystrs.no_check '[EEG LASTCOM] = pop_firpm(EEG);' catchstrs.new_and_hist];
    comfirma = [trystrs.no_check '[EEG LASTCOM] = pop_firma(EEG);' catchstrs.new_and_hist];

    % create menus if necessary
    % -------------------------
    uimenu( menu, 'Label', 'Basic FIR filter (new, default)', 'CallBack', comfirfiltnew, 'Separator', 'on', 'position', 1);
    uimenu( menu, 'Label', 'Windowed sinc FIR filter', 'CallBack', comfirws, 'position', 2);
    uimenu( menu, 'Label', 'Parks-McClellan (equiripple) FIR filter', 'CallBack', comfirpm, 'position', 3);
    uimenu( menu, 'Label', 'Moving average FIR filter', 'CallBack', comfirma, 'position', 4);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% *********************************************************************
%    Workflow for reconstruction involving variant refinement
% *********************************************************************
home; close all; clear variables;
currentFolder;
startup_mtex;
setMTEXpref('xAxisDirection','east');
setMTEXpref('zAxisDirection','outOfPlane');
setMTEXpref('FontSize',14)
ebsd = mtexdata('martensite');
%% Compute, filter and smooth grains
[grains,ebsd.grainId] = calcGrains(ebsd('indexed'),'angle',2*degree);
ebsd(grains(grains.grainSize < 3)) = [];
[grains,ebsd.grainId] = calcGrains(ebsd('indexed'),'angle',2*degree);
grains = smooth(grains,5);
%% Define and refine parent-to-child orientation relationship
job_temp = parentGrainReconstructor(ebsd,grains);
job_temp.p2c = orientation.GreningerTrojano(job_temp.csParent, job_temp.csChild);
job_temp.calcParent2Child;
%% Reconstruct parent microstructure and plot the results
job_temp.calcHyperGraph3('threshold',5*degree,'c2c','mergeSimilar','mergethreshold',8*degree,'keepGraph');
job_temp.clusterHyperGraph3('numIter',20,'inflationPower',1,'merged','mergethreshold',8*degree,'keepGraph');
job_temp.calcParentFromVote;
figure;
plot(job_temp.parentGrains,job_temp.parentGrains.meanOrientation,'linewidth',2);
%% Calculate variants and reconstruct new grains based on variants
job_temp.calcVariants;
ebsdC = job_temp.ebsdPrior(job_temp.grainsPrior(job_temp.isTransformed));
oriP = job_temp.grains(job_temp.mergeId(ebsdC.grainId)).meanOrientation;
[varIds,packIds] = calcVariantId(oriP,ebsdC.orientations,job_temp.p2c,...
                                 'variantMap', job_temp.variantMap);
varPids = [varIds,job_temp.grains(job_temp.mergeId(ebsdC.grainId)).id];
[laths,ebsdC.grainId] = calcGrains(ebsdC,'variants',varPids);
laths(ebsdC.grainId).prop.packetId = packIds;
%% Redo the reconstruction with refined child grains
job = parentGrainReconstructor(ebsdC,laths);
job.p2c = orientation.GreningerTrojano(job_temp.csParent, job_temp.csChild);
job.calcParent2Child;
%% Reconstruct parent microstructure and plot the results
job.calcHyperGraph3('threshold',5*degree,'c2c','mergeSimilar','mergethreshold',8*degree,'keepGraph');
job.clusterHyperGraph3('numIter',20,'inflationPower',1,'merged','mergethreshold',8*degree,'keepGraph');
job.calcParentFromVote;
figure;
plot(job.parentGrains,job.parentGrains.meanOrientation,'linewidth',2);
%% Clean reconstructed grains
job.mergeSimilar('threshold',7.5*degree);
job.mergeInclusions('maxSize',150);
figure;
plot(job.parentGrains,job.parentGrains.meanOrientation,'linewidth',2)
%% Variant analysis
job_temp.calcVariants;
plotMap_variants(job,'linewidth',3);
plotMap_packets(job,'linewidth',3);

author: mtex-toolbox

repo: mtex

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


text:
% Author: Paul Kienzle <pkienzle@users.sf.net> (2007)
% This program is granted to the public domain.

% -*- texinfo -*-
% @deftypefn  {Function File} {@var{y} =} downsample (@var{x}, @var{n})
% @deftypefnx {Function File} {@var{y} =} downsample (@var{x}, @var{n}, @var{offset})
% Downsample the signal, selecting every @var{n}th element.  If @var{x}
% is a matrix, downsample every column.
%
% For most signals you will want to use @code{decimate} instead since
% it prefilters the high frequency components of the signal and
% avoids aliasing effects.
%
% If @var{offset} is defined, select every @var{n}th element starting at
% sample @var{offset}.
% @seealso{decimate, interp, resample, upfirdn, upsample}
% @end deftypefn

function y = downsample(x, n, phase)

if nargin<2 || nargin>3, print_usage; end

if nargin<3
  phase = 0;
end

if phase > n - 1
  warning('This is incompatible with Matlab (phase = 0:n-1). See octave-forge signal package release notes for details.')
end

if isvector(x)
  y = x(phase + 1:n:end);
else
  y = x(phase + 1:n:end,:);
end

%!assert(downsample([1,2,3,4,5],2),[1,3,5]);
%!assert(downsample([1;2;3;4;5],2),[1;3;5]);
%!assert(downsample([1,2;3,4;5,6;7,8;9,10],2),[1,2;5,6;9,10]);
%!assert(downsample([1,2,3,4,5],2,1),[2,4]);
%!assert(downsample([1,2;3,4;5,6;7,8;9,10],2,1),[3,4;7,8]);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%DEMOHISTOGRAMRGB  One-line description here, please.
%
%   output = demoHistogramRGB(input)
%
%   Example
%   demoHistogramRGB
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2018-10-08,    using Matlab 9.4.0.813654 (R2018a)
% Copyright 2018 INRA - Cepia Software Platform.

%% Read Color image

img = Image.read('peppers.png');
figure; show(img);
print(gcf, 'peppers-show.png', '-dpng');

%% Compute histogram

figure; histogram(img);
print(gcf, 'peppers-histogram.png', '-dpng');

%% Compute line profile

lineProfile(img, [150 50], [150 300]);
print(gcf, 'peppers-lineProfile.png', '-dpng');

% overlay line profile on original image
img(149:151, 50:300, :) = 255;
write(img, 'peppers-lineOverlay.png');

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function state = AssignRandomInitialStatematrix(x,Q);

state=floor(1+Q*rand(size(x,1),size(x,1)));

repo: mcbench-benchmarks

path: github-repos/MATLAB/Sable-mcbench-benchmarks/mcbench-benchmarks-ba13b2f0296ef49491b95e3f984c7c41fccdb6d8/34784-monte-carlo-simulation-of-two-dimensional-grain-growth-code-version-no-1-basic/Monte Carlo Simulation Q-state Potts model 2D square-lattice - microstructure/AssignRandomInitialStateMatrixQPOTTS.m


text:
function [rhsQ] = CurvedEulerRHS2D(Q, time, SolutionBC, fluxtype)

% function [rhsQ] = CurvedEulerRHS2D(Q, time, SolutionBC, fluxtype)
% Purpose: compute right hand side residual for the compressible Euler
%          gas dynamics equations

Globals2D;

% 1.1 Interpolate solution to cubature nodes
cQ = zeros(cub.Ncub, K, 4);
for n=1:4, cQ(:,:,n) = cub.V*Q(:,:,n); end;

% 1.2 Evaluate flux function at cubature nodes
gamma = 1.4;
[F,G,rho,u,v,p] = EulerFluxes2D(cQ, gamma);

% 1.3 Compute volume terms (dphidx, F) + (dphidy, G)
rhsQ = zeros(Np, K, 4);
for n=1:4
  ddr = (cub.Dr')*(cub.W.*(cub.rx.*F(:,:,n) + cub.ry.*G(:,:,n)));
  dds = (cub.Ds')*(cub.W.*(cub.sx.*F(:,:,n) + cub.sy.*G(:,:,n)));
  rhsQ(:,:,n) = ddr + dds;
end

% 2.1 SURFACE TERMS (using Gauss nodes on element faces)
nx = gauss.nx; ny = gauss.ny;
mapW = gauss.mapW; mapI = gauss.mapI; mapO = gauss.mapO; mapB = gauss.mapB; mapC = gauss.mapC;

% 2.2 Interpolate solution to Gauss surface nodes
gQM = zeros(gauss.NGauss*Nfaces, K, 4); gQP = zeros(gauss.NGauss*Nfaces, K, 4);
for n=1:4
  gQ = gauss.interp*Q(:,:,n);
  gQM(:,:,n) = gQ(gauss.mapM);  gQP(:,:,n) = gQ(gauss.mapP);
end

% 2.3 Apply boundary conditions to '+' traces
if(~isempty(SolutionBC))
  gQP = feval(SolutionBC, gauss.x, gauss.y, gauss.nx, gauss.ny, gauss.mapI, ...
                 gauss.mapO, gauss.mapW, gauss.mapC, gQP, time);
end

% 2.4 Evaluate surface flux functions with stabilization
switch fluxtype
  case {'LF'}
    [flux] = EulerLF2D (nx, ny, gQM, gQP, gamma);
  case {'Roe'}
    [flux] = EulerRoe2D(nx, ny, gQM, gQP, gamma);
  case {'HLL'}
    [flux] = EulerHLL2D(nx, ny, gQM, gQP, gamma);
  case {'HLLC'}
    [flux] = EulerHLLC2D(nx, ny, gQM, gQP, gamma);
end

% 2.5 Compute surface integral terms
for n=1:4
  rhsQ(:,:,n) = rhsQ(:,:,n) - gauss.interp'*(gauss.W.*flux(:,:,n));
end

% 3.1 Multiply by inverse mass matrix
for n=1:4
  % 3.1.a Multiply straight sided elements by inverse mass matrix
  rhsQ(:,straight, n) = V*V'*(rhsQ(:,straight, n)./J(:,straight));

  % 3.1.b Multiply curvilinear elements by custom inverse mass matrices
  Ncurved = length(curved);
  for m=1:Ncurved
    k = curved(m);
    mmCHOL = cub.mmCHOL(:,:,k);
    rhsQ(:,k,n) = mmCHOL\(mmCHOL'\rhsQ(:,k,n));
  end
end
return;

repo: Aero-matlab