text:
function a_inverse = r8ut_inverse ( n, a )

%*****************************************************************************80
%
%% R8UT_INVERSE computes the inverse of a R8UT matrix.
%
%  Discussion:
%
%    The R8UT storage format is used for an M by N upper triangular
%    matrix.  The format stores all M*N entries, even those which are zero.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    A Nijenhuis and H Wilf,
%    Combinatorial Algorithms,
%    Academic Press, 1978, second edition,
%    ISBN 0-12-519260-6
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real A(N,N), the R8UT matrix to be inverted.
%
%    Output, real A_INVERSE(N,N), the inverse of the upper triangular matrix.
%

%
%  Check.
%
  for i = 1 : n
    if ( a(i,i) == 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R8UT_INVERSE - Fatal error!\n' );
      fprintf ( 1, '  Zero diagonal element.\n' );
      error ( 'R8UT_INVERSE - Fatal error!' );
    end
  end

  a_inverse(1:n,1:n) = a(1:n,1:n);

  for j = n : -1 : 1

    for i = n : -1 : 1

      if ( j < i )

        a_inverse(i,j) = 0.0;

      elseif ( i == j )

        a_inverse(i,j) = 1.0 / a_inverse(i,j);

      elseif ( i < j )

        a_inverse(i,j) = - ( a_inverse(i,i+1:j) * a_inverse(i+1:j,j) ) ...
          / a_inverse(i,i);

      end

    end
  end

  return
end

repo: jburkardt-m


text:
clear;
clc;

s = iio_sys_obj_matlab; % Constructor
s.ip_address = '10.66.99.200';
s.dev_name = 'adrv9009';
s.in_ch_no = 2;
s.out_ch_no = 2;
s.in_ch_size = 8192;
s.out_ch_size = 8192;

s = s.setupImpl();

input = cell(1, s.in_ch_no + length(s.iio_dev_cfg.cfg_ch));
Fs = 245.76e6;
Fc = 1e6;
Bw = 200e6;
t = 1/Fs:1/Fs:s.in_ch_size/Fs;
for i=1:s.in_ch_no
    input{i} = sin(2*pi*Fc*t+(i-1)*pi/2)*2^15;
end

input{s.getInChannel('TRX_LO_FREQ')} = 2.45e9;
input{s.getInChannel('RX1_RF_BANDWIDTH')} = Bw;
input{s.getInChannel('RX1_GAIN_MODE')} = 'slow_attack';
input{s.getInChannel('RX1_GAIN')} = 0;
input{s.getInChannel('RX2_GAIN_MODE')} = 'slow_attack';
input{s.getInChannel('RX2_RF_BANDWIDTH')} = Bw;
input{s.getInChannel('RX2_GAIN')} = 0;
input{s.getInChannel('TX_RF_BANDWIDTH')} = Bw;
input{s.getInChannel('TX_GAIN')} = -20;

for i = 1:20
output = stepImpl(s, input);
end

s.releaseImpl();

figure % new figure
ax1 = subplot(2,1,1); % top subplot
ax2 = subplot(2,1,2); % bottom subplot

plot(ax1,output{1});
title(ax1,'I');
xlabel('Sample');
ylabel('Amplitude');

plot(ax2,output{2});
title(ax2,'Q');
xlabel('Sample');
ylabel('Amplitude');

author: analogdevicesinc


text:
%PLOTO Plot objects as 1-D functions of the feature number
%
%   [HH HO HC] = PLOTO(A,N)
%
% INPUT
%   A   Dataset
%   N   Integer
%
% OUTPUT
%   HH  Lines handles
%   HO  Object identifier handles
%   HC  Class number handles
%
% DESCRIPTION
% Produces 1-D function plots for all the objects in dataset A. The plots
% are organised as subplots, N on a row. Default is the squareroot of the
% number of objects. Object identifiers and class numbers are written in
% the correspopnding plots.
%
% See also DATASETS

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function [h_out1,h_out2,h_out3] = ploto(a,p)

	if nargin < 2, p = []; end
	[m,k,c] = getsize(a);
	nlab = getnlab(a);

	% Define the color for each of the classes:
	if c == 2
		map = [0 0 1; 1 0 0];
	else
		map = hsv(c);
	end

	% Make subplots for each object, so a grid of p x q subplots is
	% defined
	h = [];
	if ~isempty(p)
		q = ceil(m/p);
	elseif m > 3
		p = ceil(sqrt(m)); q = ceil(m/p);
	else
		p = m; q = 1;
	end
	% Get the object labels
	labs = getlabels(a);
	ymin = min(a.data(:));
	ymax = max(a.data(:));
	V = [1 k ymin ymax];
	% Make the plot for each of the objects:
	h = [];
	ho = [];
	hc = [];
	s = sprintf('Plot %i objects: ',m);
	prwaitbar(m,s);
	for j = 1:m
		if isdatafile(a) | 1
			prwaitbar(m,j,[s int2str(j)]);
			b = +prdataset(a(j,:));
			ymin = min(b);
			ymax = max(b);
			k = length(b);
			V = [1 k ymin ymax];
		else
			b = +a(j,:);
		end
		% Create the subplots with the correct sizes:
		subplot(q,p,j)
		hh = plot(b);
		set(gca,'xtick',[]);
		set(gca,'ytick',[]);
		axis(gca,V);
		ho = [ho text(2,ymax-0.15*(ymax-ymin),getident(a(j,:),'string'))];
		hc = [hc text(3*k/4,ymax-0.15*(ymax-ymin),num2str(nlab(j)))];
		h = [h hh];
		hold on
	end
	prwaitbar(0);

	% The last details to take care of:
	if nargout > 0
		h_out1 = h;
		h_out2 = ho;
		h_out3 = hc;
	end

	return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function fmri_EccBias_GLMcontrasts

% EccBias scans processing code
%
% USER IMPORTANT!!!: Before running this script, make sure that all EccBias
% scans have the correct annotation fields and scan groups (are grouped
% correctly). Otherwise, fmri_EccBias_checkDataTypesForGLM will fail.
%
% 1. Set GLM parameters for all scans
% 2. Check EccBias localizer experiments for:
%    a. Datatype and scan numbers
%    b. Record parfiles
%    c. Make sure scan groupings are correct
% 3. Run GLM for EccBias experiments
% 4. Run contrasts for the EccBias experiments
%
% Calls the functions:
% * fmri_EccBias_setGLMparams (subfunction)
% * fmri_EccBias_checkDataTypesForGLM
% * fmri_EccBias_computeContrasts
%
%
% AL 09/17/08
% Modified 06/18/2008: there are calls to er_setParams and saveSessions all
% over the place, which results in the params for the MotionComp_RefScan1
% being changed for the worse. I'm not sure how to efficiently prevent this
% from happening, so instead, I let it happen, and just fix the ANNOTATION
% and PARFILE fields at the end, since accuracy here is important.
%
% Code taken from fmri_FFA_motionCompFix

clear all

% Set directory and subject list
% Alternatively: subs = {'sess1','sess2','sess3',...};
if ispc
   fmriDir = 'W:\projects\Kids\fmri\eccbias_faces_places\';
else
   fmriDir = '/biac1/kgs/projects/Kids/fmri/eccbias_faces_places/';
end
cd(fmriDir); s=dir('*0*'); subs={s.name};

% Start a log text file to document successes and failures in preprocessing
dateAndTime=datestr(now); dateAndTime(12)='_'; dateAndTime(15)='h';dateAndTime(18)='m';
logDir=fullfile(fmriDir,'logs');
if ~isdir(logDir), mkdir(logDir), end;
logFile = fullfile(fmriDir,'logs',['EccBias_Log_' dateAndTime '.txt']);
fid=fopen(logFile,'w');
[junk,logFileName]=fileparts(logFile);
scriptName = 'fmri_EccBias_GLMcontrasts.m';
startTime = clock;

% 1. Set GLM parameters
params=fmri_EccBias_setGLMparams(fid,logFileName,scriptName); % Subfunction within this script

for ii=1:length(subs)

    % Get rid of all mrVista globals in workspace. Put this at the top,
    % just in case we break from the rest of the loop.
    clear dataTYPES mrSESSION vANATOMYPATH HOMEDIR hI newDt newScan

    thisDir = fullfile(fmriDir,subs{ii}); cd(thisDir);
    fprintf(fid,'\n ------------------------------------------ \n');
    fprintf(fid,'Processing %s\n\n',thisDir);
    fprintf('Processing %s\n\n',thisDir);

    %% Check if the subject already has been processed by script (and thus
    %% should not be re-entered to this script). If they do, skip to the
    %% next subject by using CONTINUE.
    load mrSESSION.mat
    theglm=find(strcmp({dataTYPES.name},'GLMs'));
     if ~isempty(theglm)
        checkScans=length(dataTYPES(theglm).eventAnalysisParams);
        alreadyDone=[];
        for xx=1:checkScans
            tmp=(strfind(lower(dataTYPES(theglm).eventAnalysisParams(xx).annotation),'eccbias'));
            alreadyDone=[alreadyDone tmp];
        end
        if(~isempty(alreadyDone))
            fprintf(fid,'Error: Subject already processed (has GLM with ''Eccbias'' in annotation field), skipping... \n\n');
            fprintf('Error: Subject already processed (has GLM with ''Eccbias'' in annotation field), skipping... \n\n');
            continue
        end
    end

    %% 2a/b/c. Find datatype number for MotionComp_RefScan1, find scan
    %% numbers for EccBias localizer scans, record parfiles for these scans and
    %% make sure scan groupings are correct. If everything checks out,
    %% proceed with GLM. If not, SKIP THIS SUBJECT.
    [go,s1,s2]=fmri_Eccbias_checkDataTypesForGLM(fid,dataTYPES);

    % If there was a problem in step 3, skip GLM and proceed to next
    % subject
    if go==0
        continue

    %% 3. Otherwise, run GLM (to debug: fid=1)
    elseif go==1
        try
            % Initialize a hidden view of the 'inplane' data (like invoking
            % mrVista without a GUI, dataTYPES/etc in workspace)
            hI = initHiddenInplane(s1.mc, s1.scan);
            params.parfiles={s1.parfile, s2.parfile};
            er_setParams(hI, params); % write params as eventAnalysisParams fields
            tic; [hI, ecc_scan] = applyGlm(hI,s1.mc,s1.scangroup,params);
            time=toc;
            fprintf(fid,'GLM COMPLETED -- dt %d, scan %d. Time: %d m %2.2f s \n\n',...
                hI.curDataType,ecc_scan,floor(time/60),mod(time,60));
            fprintf('GLM COMPLETED -- dt %d, scan %d. Time: %d m %2.2f s \n\n',...
                hI.curDataType,ecc_scan,floor(time/60),mod(time,60));
            glmOK=1;
        catch
            fprintf(fid,'FAILURE TO SUCCESSFULLY RUN GLM \n');
            fprintf('FAILURE TO SUCCESSFULLY RUN GLM \n');
            theerror=lasterror;
            fprintf(fid,'%s \n\n',theerror.message);
            fprintf('%s \n\n',theerror.message);
            glmOK=0;
        end

        %% 4. Compute contrast maps
        if glmOK

            try
                fmri_EccBias_computeContrasts(hI, ecc_scan)
                fprintf(fid,'Computed contrasts using fmri_EccBias_computeContrasts \n');
                fprintf('Computed contrasts using fmri_EccBias_computeContrasts \n');
            catch
                fprintf(fid,'FAILURE TO SUCCESSFULLY COMPUTE CONTRASTS \n');
                fprintf('FAILURE TO SUCCESSFULLY COMPUTE CONTRASTS \n');
                theerror=lasterror;
                fprintf(fid,'%s \n\n',theerror.message);
                fprintf('%s \n\n',theerror.message);

            end
        end
    end

    % Fix the annotation and parfile fields for the MotionComp_RefScan
    % dataTYPE, since these are undesirably changed (collateral damage).
    dataTYPES(s1.mc).eventAnalysisParams(s1.scan).annotation='EccBias_run1';
    dataTYPES(s2.mc).eventAnalysisParams(s2.scan).annotation='EccBias_run2';

    par1=dataTYPES(s1.mc).scanParams(s1.scan).parfile;
    par2=dataTYPES(s2.mc).scanParams(s2.scan).parfile;

    dataTYPES(s1.mc).eventAnalysisParams(s1.scan).parfiles=par1;
    dataTYPES(s2.mc).eventAnalysisParams(s2.scan).parfiles=par2;
    save(fullfile(thisDir,'mrSESSION.mat'),'mrSESSION','dataTYPES', '-append');
    clear par1 par2

end

totalTime=etime(clock,startTime);

fprintf(fid,'\n ------------------------------------------ \n');
fprintf(fid,'Total running time for script: %f minutes \n',totalTime/60);
fprintf('Total running time for script: %f minutes \n',totalTime/60);

fclose(fid); % close out the log file

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function params=fmri_EccBias_setGLMparams(fid,logFileName,scriptName)

% Set GLM parameters: We set the PARAMS struct fields to the values on the
% wiki, which are further explained there. Also we print the values to the
% log. http://vpnl.stanford.edu/internal/wiki/index.php/ROIanalysis#GLMs

fprintf(fid,'\n ------------------------------------------ \n');
fprintf(fid,'\n GLM PARAMETERS \n');
fprintf(fid,' ------------------------------------------ \n');

params=er_defaultParams % Initialize param struct with defaults
p.al='params.alpha=0.0500'; eval(p.al); fprintf(fid,'%s \n',p.al);
p.at='params.ampType=''betas'''; eval(p.at); fprintf(fid,'%s \n',p.at);
p.an='params.annotation=''Eccbias Runs 1&2'''; eval(p.an); fprintf(fid,'%s \n',p.an);
p.ap='params.assignParfiles=0'; eval(p.ap); fprintf(fid,'%s \n',p.ap);
p.bp='params.bslPeriod=[-6 -5 -4 -3 -2 -1 0]'; eval(p.bp); fprintf(fid,'%s \n',p.bp);
p.dt='params.detrend=1'; eval(p.dt); fprintf(fid,'%s \n',p.dt);
p.df='params.detrendFrames=20;'; eval(p.df); fprintf(fid,'%s \n',p.df);
p.ea='params.eventAnalysis=1'; eval(p.ea); fprintf(fid,'%s \n',p.ea);
p.eb='params.eventsPerBlock=6'; eval(p.eb); fprintf(fid,'%s \n',p.eb);
p.fp='params.framePeriod=2'; eval(p.fp); fprintf(fid,'%s \n',p.fp);
p.gt='params.glmHRF=3'; eval(p.gt); fprintf(fid,'%s \n',p.gt);
p.gp='params.glmHRF_params=[6 16 1 1 6 0 24]'; eval(p.gp); fprintf(fid,'%s \n',p.gp);
p.gw='params.glmWhiten=0'; eval(p.gw); fprintf(fid,'%s \n',p.gw); % temporal autocorrelation
p.ic='params.inhomoCorrect=1'; eval(p.ic); fprintf(fid,'%s \n',p.ic);
p.lp='params.lowPassFilter=0'; eval(p.lp); fprintf(fid,'%s \n',p.lp);
p.nb='params.normBsl=1'; eval(p.nb); fprintf(fid,'%s \n',p.nb);
p.od='params.onsetDelta=0'; eval(p.od); fprintf(fid,'%s \n',p.od);
p.pf='params.parfiles='''''; eval(p.pf); fprintf(fid,'%s \n',p.pf); % fill in after checkDataTypes...
p.pp='params.peakPeriod=[4 5 6 7 8 9 10 11 12 13 14]'; eval(p.pp); fprintf(fid,'%s \n',p.pp);
p.sh='params.setHRFParams=0'; eval(p.sh); fprintf(fid,'%s \n',p.sh);
p.sc='params.snrConds=1'; eval(p.sc); fprintf(fid,'%s \n',p.sc);
p.tn='params.temporalNormalization=0'; eval(p.tn); fprintf(fid,'%s \n',p.tn);
p.tw='params.timeWindow=[-6:1:24]'; eval(p.tw); fprintf(fid,'%s \n',p.tw);
p.log=['params.logFile=''' logFileName '''']; eval(p.log); fprintf(fid,'%s \n',p.log)
p.script=['params.scriptName=''' scriptName '''']; eval(p.script); fprintf(fid,'%s \n',p.script)

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%kgtubeness 'Replace Lines with Tesselated Cylinders'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros gtubeness.pane file
%
% Parameters:
% InputFile: i 'Input Geometry ', required: 'Input geometry object.'
% Double: scale 'Radius Scale', default: 0.1: 'A scalar which indicates the radius of the tube.'
% Integer: subdivisions 'Subdivisions', default: 3: 'The number of sides for each tube will be 2^(subdivisions - 1).'
% OutputFile: o 'Output Geometry', required: 'Tube-ified geometry object'
%
% Example: o = kgtubeness(i, {'i','';'scale',0.1;'subdivisions',3;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% gtubeness - Replace Lines with Tesselated Cylinders
%
%  DESCRIPTION
% Filter for geometry objects in which line segments will be replaced with
% tubes.  Non-line-segment primitives will be passed through unchanged.
%
% Currently the only "cap_style" options that work are "None" and "Round caps".
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1996,1997 , The Regents of the University of California.  All rights reserved.
%

function varargout = kgtubeness(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,..] = kgtubeness(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'scale', 0.1;'subdivisions', 3;'o', '__output'};
maxval={0,1,6,0};
minval={0,1,2,0};
istoggle=[0,1,1,0];
was_set=istoggle * 0;
paramtype={'InputFile','Double','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 'gtubeness"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [ homos ] = NewWarping( pa_, pb_, H, W, qH, qW, lambda)
% A warpped up wersion of as-similar-as-possible warping, with pre-warp.
% set PREWARP to false if the output is badly distorted.
    PREWARP = true;
    nP = length(pa_);
    if length(pb_) ~= nP
        error('Points Numbers met Matching!');
    end
    if PREWARP
        [preH, ~] = ransacfithomography(pa_', pb_', 0.001);
        whilecount = 0;
        while isnan(preH(1, 1)) && whilecount < 100
            [preH, ~] = ransacfithomography(pa_', pb_', 0.001);
            whilecount = whilecount + 1;
        end
        if whilecount == 100
            error('?') ;
        end
    else
        preH = eye(3);
        pa = pa_;
    end
    pbWarp = preH \ [pb_' ; ones(1, nP)];
    pbWarp(1, :) = pbWarp(1, :) ./ pbWarp(3, :);
    pbWarp(2, :) = pbWarp(2, :) ./ pbWarp(3, :);
    pbWarp = pbWarp(1:2, :)';

    diff = sum((pa_ - pbWarp) .* (pa_ - pbWarp), 2) < 1000;
    valid = pbWarp(:, 1) > 0 & pbWarp(:, 1) < W & pbWarp(:, 2) > 0 & pbWarp(:, 2) < H;
    pa = pa_(valid & diff, :);
    pbWarp = pbWarp(valid & diff, :);

    asap = AsSimilarAsPossibleWarping(H, W, qW, qH, lambda);
    asap.SetControlPts(pa, pbWarp);
    asap.Solve();

% -----DEBUG-SCRIPT------
% use it to check the warping result
%     e = asap.CalcError();

%     grid = asap.Warp(ones(H, W, 3) * 255, 400);
%     imshow(grid);

%     if e > 10
%         disp('?') ;
%     end

    homos2 = asap.CalcHomos();
    homos = homos2;
    for row = 1:H/qH
        for col = 1:W/qW
            tempH = preH * squeeze(homos2(row, col, :, :));
            homos(row, col, :, :) = tempH ./ tempH(3, 3);
        end
    end
end

author: SuTanTank

repo: VideoStitchingViaShakinessRemoving


text:
%%  mnist data
clear all; close all; clc;
load mnist_uint8;
x = cell(100, 1);
N = 600;
for i = 1 : 100
    x{i}{1} = reshape(train_x(((i - 1) * N + 1) : (i) * N, :), N, 28, 28) * 255;
end
%% ex 1
scae = {
    struct('outputmaps', 10, 'inputkernel', [1 5 5], 'outputkernel', [1 5 5], 'scale', [1 2 2], 'sigma', 0.1, 'momentum', 0.9, 'noise', 0)
};

opts.rounds     = 1000;
opts.batchsize  =    1;
opts.alpha      = 0.01;
opts.ddinterval =   10;
opts.ddhist     =  0.5;
scae = scaesetup(scae, x, opts);
scae = scaetrain(scae, x, opts);
cae = scae{1};

%Visualize the average reconstruction error
plot(cae.rL);

%Visualize the output kernels
ff=[];
for i=1:numel(cae.ok{1});
    mm = cae.ok{1}{i}(1,:,:);
    ff(i,:) = mm(:);
end;
figure;visualize(ff')

author: rasmusbergpalm

repo: DeepLearnToolbox


text:
function f = p37_f ( n, x )

%*****************************************************************************80
%
%% P37_F evaluates the objective function for problem 37.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Zbigniew Michalewicz,
%    Genetic Algorithms + Data Structures = Evolution Programs,
%    Third Edition,
%    Springer Verlag, 1996,
%    ISBN: 3-540-60676-9,
%    LC: QA76.618.M53.
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the argument of the objective function.
%
%    Output, real F, the value of the objective function.
%
  arg = - ( x(1) - pi )^2 - ( x(2) - pi )^2;
  f = - cos ( x(1) ) * cos ( x(2) ) * exp ( arg );

  return
end

author: johannesgerer


text:
%  plots prior and, if information is available, posterior densities
%
%  ::
%
%    [h,legend_]=prior_posterior(ss,varargin)
%
%  Args:
%
%     - **ss** [struct]: structure containing the relevant information to plot
%       and more specifically
%
%       - **x_kdens** [vector]: x-axis values for the posterior density
%       - **f_kdens** [vector]: y-axis values for the posterior density
%       - **x_prior** [vector]: x-axis values for the prior density
%       - **f_prior** [vector]: y-axis values for the prior density
%       - **mean_sim** [scalar]: mean of the posterior simulation
%       - **post_mode** [scalar]: value at the posterior maximization mode
%       - **post_mode_sim** [scalar]: value at the posterior simulation mode
%       - **tex_name** [char]: name of the parameter
%
%     - **varargin** [pairwise arguments]: standard plotting arguments for
%       matlab
%
%  Returns:
%     :
%
%     - **h** [handle]: handle for the plot
%     - **legend_** [cellstr]: names of the lines in the plot
%     - **tex_name** [char]: name of the parameter
%
%  Note:
%
%     - Only the prior density is plotted if no posterior information is
%       available.
%
%

author: jmaih


text:
% Measure PSS correlations in the presence of a frequency offset.
%n_id_2=0;

% Copyright 2012 Evrytania LLC (http://www.evrytania.com)
%
% Written by James Peroulas <james@evrytania.com>
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU Affero 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 Affero General Public License for more details.
%
% You should have received a copy of the GNU Affero General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

freq_off_set=linspace(-30e3,30e3,1001);
time_off_set=-3000:3000;

n_freqs=length(freq_off_set);
n_times=length(time_off_set);

pss_td=NaN(3,144+2048+0*2048);
for n_id_2=0:2
  pss_freq=pss(n_id_2);
  pss_freq=[0 pss_freq(32:end) zeros(1,2048-73+5+5) pss_freq(1:31)];
  pss_td(n_id_2+1,145:end)=[idft(pss_freq)*sqrt(2048/62) zeros(1,0*2048)];
  pss_td(n_id_2+1,1:144)=pss_td(n_id_2+1,end-143:end);
  pss_td(n_id_2+1,:)=pss_td(n_id_2+1,:)/sqrt(sigpower(pss_td(n_id_2+1,:)));
end

% Frequency domain correlations
figure(1);
log_xc_freq=NaN(3,n_freqs);
for n_id_2=0:2
  for t=1:n_freqs
    freq=freq_off_set(t);
    log_xc_freq(n_id_2+1,t)=sum(conj(pss_td(n_id_2+1,:)).*fshift(pss_td(n_id_2+1,:),freq/(fs_lte/2)))/(2048+144);
  end
end
plot(freq_off_set,db20(abs(transpose(log_xc_freq))));
%ylim([-20 5]);
drawnow;
zgo;

% Time domain auto correlations
pss_td_ext=[pss_td zeros(3,(2048+144)*2)];
log_xc_td=NaN(3,n_times);
figure(2);
for n_id_2=0:2
  for t=1:n_times
    to=time_off_set(t);
    log_xc_td(n_id_2+1,t)=sum(conj(pss_td_ext(n_id_2+1,:)).*tshift(pss_td_ext(n_id_2+1,:),to))/(2048+144);
  end
end
plot(time_off_set,db20(abs(transpose(log_xc_td))));
%ylim([-50 5]);
drawnow;
zgo;

% Time domain cross correlations
log_xc_td_cross=NaN(3,n_times);
log_xc_td_cross_max16=NaN(3,n_times);
figure(3);
xc_set=[1 2;1 3;2 3]-1;
for k=1:3
  n_id_2_1=xc_set(k,1);
  n_id_2_2=xc_set(k,2);
  for t=1:n_times
    to=time_off_set(t);
    log_xc_td_cross(k,t)=sum(conj(pss_td_ext(n_id_2_1+1,:)).*tshift(pss_td_ext(n_id_2_2+1,:),to))/(2048+144);
  end
  for t=1:n_times
    log_xc_td_cross_max16(k,t)=max(log_xc_td_cross(k,max([t-32 1]):min([t+80 n_times])));
  end
end
plot(time_off_set,db20(abs(transpose(log_xc_td_cross))), ...
time_off_set,db20(abs(transpose(log_xc_td_cross_max16))));
legend('0 1','0 2','2 3','location','se');
%ylim([-50 5]);
drawnow;
zgo;

author: JiaoXianjun

repo: rtl-sdr-LTE


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

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

% Last Modified by GUIDE v2.5 03-Jan-2014 18:04:48

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

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

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

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

    hRtsMainGUI = varargin{1};
    vesselGUID = varargin{2};
    handles.mainGUI = hRtsMainGUI;
    handles.vesselGUID = vesselGUID;

    if(ishandle(handles.resSumLabel1))
        rts_initRTSPropPanel(hRtsMainGUI, hObject, handles.clockLabel, handles.connStatusLabel, vesselGUID, handles.resUITable1, handles.resSumLabel1, ...
                             handles.resUITable2, handles.resSumLabel2, handles.resourceCombo1, handles.resourceCombo2, handles.mainEngineUITable,...
                             handles.rcsEngineUITable, handles.resourcePanel1, handles.resourcePanel2);
    end

    % Update handles structure
    guidata(hObject, handles);

    % UIWAIT makes rts_PropConsole wait for user response (see UIRESUME)
    % uiwait(handles.rts_PropConsole);
end

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

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

% --- Executes on selection change in resourceCombo2.
function resourceCombo2_Callback(hObject, eventdata, handles)
% hObject    handle to resourceCombo2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns resourceCombo2 contents as cell array
%        contents{get(hObject,'Value')} returns selected item from resourceCombo2
end

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

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

% --- Executes on selection change in resourceCombo1.
function resourceCombo1_Callback(hObject, eventdata, handles)
% hObject    handle to resourceCombo1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns resourceCombo1 contents as cell array
%        contents{get(hObject,'Value')} returns selected item from resourceCombo1
end

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

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

% --- Executes on selection change in resourceCombo2.
function popupmenu4_Callback(hObject, eventdata, handles)
% hObject    handle to resourceCombo2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns resourceCombo2 contents as cell array
%        contents{get(hObject,'Value')} returns selected item from resourceCombo2
end

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

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

% --- Executes when user attempts to close rts_PropConsole.
function rts_PropConsole_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to rts_PropConsole (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: delete(hObject) closes the figure
    rts_pageCloseStopTimers(hObject);
    delete(hObject);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%%*****************************************************************************
%% HSDsqlpmisc:
%% produce infeasibility certificates if appropriate
%%
%% Input: X,y,Z are the original variables, not the HSD variables.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004.
%%*****************************************************************************

   function [X,y,Z,resid,reldist,param,msg] = HSDsqlpmisc(blk,At,C,b,X,y,Z,permZ,param);

   obj         = param.obj;
   relgap      = param.relgap;
   prim_infeas = param.prim_infeas;
   dual_infeas = param.dual_infeas;
   ZpATynorm   = param.ZpATynorm;
   inftol      = param.inftol;
   m0          = param.m0;
   indeprows   = param.indeprows;
   termcode    = param.termcode;
   AX          = param.AX;
   normX0      = param.normX0;
   normZ0      = param.normZ0;
   printlevel  = param.printlevel;
%%
   resid = []; reldist = []; msg = [];
   Anorm = ops(At,'norm'); xnorm = ops(X,'norm'); ynorm = norm(y);
%%
   if (termcode <= 0)
      %%
      %% To detect near-infeasibility when the algorithm provides
      %% a "better" certificate of infeasibility than of optimality.
      %%
      err = max([prim_infeas,dual_infeas,relgap]);
      iflag = 0;
      if (obj(2) > 0)
         homRd = ZpATynorm/obj(2);
         if (homRd < 1e-2*sqrt(err*inftol))
            iflag = 1;
            termcode = 1;
            param.termcode = 1;
         end
      elseif (obj(1) < 0)
         homrp = norm(AX)/(-obj(1));
         if (homrp < 1e-2*sqrt(err*inftol))
            iflag = 1;
            termcode = 2;
            param.termcode = 2;
         end
      end
   end
   if (termcode == 1)
      rby = 1/(b'*y); y = rby*y; Z = ops(Z,'*',rby);
      resid = ZpATynorm * rby;
      reldist = ZpATynorm/(Anorm*ynorm);
      msg = 'Stop: primal problem is suspected of being infeasible';
      if (printlevel); fprintf('\n  %s',msg); end
   end
   if (termcode == 2)
      tCX = blktrace(blk,C,X);
      X = ops(X,'*',1/(-tCX));
      resid = norm(AX) /(-tCX);
      reldist = norm(AX)/(Anorm*xnorm);
      msg = 'Stop: dual problem is suspected of being infeasible';
      if (printlevel); fprintf('\n  %s',msg); end
   end
   if (termcode == 3)
      maxblowup = max(ops(X,'norm')/normX0,ops(Z,'norm')/normZ0);
      msg = sprintf('Stop: primal or dual is diverging, %3.1e',maxblowup);
      if (printlevel); fprintf('\n  %s',msg); end
   end
   [X,Z] = unperm(blk,permZ,X,Z);
   if ~isempty(indeprows)
      ytmp = zeros(m0,1);
      ytmp(indeprows) = y;
      y = ytmp;
   end
%%*****************************************************************************
%% unperm: undo the permutations applied in validate.
%%
%% [X,Z,Xiter,Ziter] = unperm(blk,permZ,X,Z,Xiter,Ziter);
%%
%% undoes the permutation introduced in validate.
%% can also be called if Xiter and Ziter have not been set as
%%
%% [X,Z] = unperm(blk,permZ,X,Z);
%%
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 2 Feb 01
%%*****************************************************************************

  function [X,Z] = unperm(blk,permZ,X,Z);
%%
  for p = 1:size(blk,1)
     if (strcmp(blk{p,1},'s') & ~isempty(permZ{p}))
        per = permZ{p};
        X{p} = X{p}(per,per);
        Z{p} = Z{p}(per,per);
     end
  end
%%*****************************************************************************

repo: REST


text:
function [B,m,ps,Population] = UpdateParameter(Problem,Population)
% Update the parameters in ARSBX

%--------------------------------------------------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB Platform
% for Evolutionary Multi-Objective Optimization [Educational Forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    PopDec = Population.decs;
    Flag   = Population.adds(zeros(length(Population),1));
    Ori    = sum(Flag == 1);
    Eig    = sum(Flag == 2);
    ps     = 1/(1+exp(-Problem.M*sqrt(Problem.D)*((Ori+1)/(Eig+Ori+2)-0.5)*Problem.FE/Problem.maxFE));
    C        = cov(PopDec);
    [B,E]    = eig(C);
    E        = diag(E);
    E        = sqrt(E);
    [~,Rank] = sort(E,'descend');
    B        = B(:,Rank);
    m        = mean(PopDec);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


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

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

Vfrac_initial = 1;
optimality_initial = 1e-2;
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 denoise_PRL_2010_main
%
% This is a demo program of the paper J. Tian, W. Yu, and L. Ma, "AntShrink: Ant
% colony optimization for image shrinkage," Pattern Recognition Letters,
% Vol. 31, Oct. 2010, pp. 1751-1758.
%
% Contact: eejtian@gmail.com
%
% Note that the PSNR values could be slightly different with that reported
% in paper due to the random number generator used to generate noisy image.
%
% The input image should have a square size.
%
% Acknowledgement: This program needs Wavelet transform toolbox, which is
% downloaded from http://www-stat.stanford.edu/~wavelab/.

clear all; close all; clc;

% Load the ground truth image
img_truth = double(imread('barbara_truth.bmp'));
[nRow, nColumn] = size(img_truth);

% Use the ground truth image to generate the noisy image
noise_sig_truth = 10; % sigma_n used in the paper. This parameter is adjusted by the user.
noise_mu = 0;
img_noisy = img_truth + randn(size(img_truth)) .* noise_sig_truth + noise_mu;

% wavelet parameters
wbase = 'Daubechies';
mom = 8;
dwt_level = 5; %note that here, dwt_scale means the decomposition level of the DWT
[n,J] = func_quadlength(img_truth);
L = J-dwt_level;%here, L means the size of the coarsest level, 2^L

win_size = 2;
img_denoised = zeros(size(img_noisy));

% Since this is time consuming approach, divide the image into four parts, then
% conduct denoising for each part
for ii=1:4
    win_size=2;

    switch ii
        case 1
            img_denoised(1:end/2,1:end/2) = func_ACOShrink(img_noisy(1:end/2,1:end/2), wbase, mom, dwt_level, win_size);
        case 2
            img_denoised(end/2+1:end,1:end/2) = func_ACOShrink(img_noisy(end/2+1:end,1:end/2), wbase, mom, dwt_level, win_size);
        case 3
            img_denoised(1:end/2,end/2+1:end) = func_ACOShrink(img_noisy(1:end/2,end/2+1:end), wbase, mom, dwt_level, win_size);
        case 4
            img_denoised(end/2+1:end,end/2+1:end) = func_ACOShrink(img_noisy(end/2+1:end,end/2+1:end), wbase, mom, dwt_level, win_size);
    end

end

% Calculate the PSNR performance
fprintf('PSNR=%.2fdB\n', func_psnr_gray(img_truth, img_denoised));

% Write the output image
imwrite(uint8(img_denoised), 'barbara_denoised.bmp','bmp');

%-------------------------------------------------------------------------
%------------------------------Inner Function ----------------------------
%-------------------------------------------------------------------------
% Main algorithm of proposed AntShrink algorithm
function x_out= func_ACOShrink(x_in, wbase, mom, dwt_level, win_size)

[nrow, ncol] = size(x_in);
L = log2(size(x_in,2))-dwt_level;

% Estimate the noise_sigma from the noisy signal
qmf = func_MakeONFilter(wbase,mom);
[temp, coef] = func_NormNoise_2d(x_in, qmf);
noise_sigma = 1/coef;

wx  = func_FWT2_PO(x_in, L, qmf);
[n,J] = func_dyadlength(wx);
ws = wx;

rr = (meshgrid(1:nrow))'; cc = meshgrid(1:ncol);

ant_search_range_row_min = rr;
ant_search_range_row_max = rr;
ant_search_range_col_min = cc;
ant_search_range_col_max = cc;

for j=(J-1):-1:L
    [t1,t2] = func_dyad2HH(j);
    ant_search_range_row_min(t1,t2) = min2(t1);
    ant_search_range_row_max(t1,t2) = max2(t1);
    ant_search_range_col_min(t1,t2) = min2(t2);
    ant_search_range_col_max(t1,t2) = max2(t2);
    [t1,t2] = func_dyad2HL(j);
    ant_search_range_row_min(t1,t2) = min2(t1);
    ant_search_range_row_max(t1,t2) = max2(t1);
    ant_search_range_col_min(t1,t2) = min2(t2);
    ant_search_range_col_max(t1,t2) = max2(t2);
    [t1,t2] = func_dyad2LH(j);
    ant_search_range_row_min(t1,t2) = min2(t1);
    ant_search_range_row_max(t1,t2) = max2(t1);
    ant_search_range_col_min(t1,t2) = min2(t2);
    ant_search_range_col_max(t1,t2) = max2(t2);
end

data_var = func_signal_variance_estimation(wx, noise_sigma, win_size,ant_search_range_row_min,ant_search_range_row_max,ant_search_range_col_min,ant_search_range_col_max,J,L);

for j=(J-1):-1:L
    [t1,t2] = func_dyad2HH(j);
    ws(t1,t2) = wx(t1,t2) .* data_var(t1,t2) ./ (data_var(t1,t2) + noise_sigma.^2);
    [t1,t2] = func_dyad2HL(j);
    ws(t1,t2) = wx(t1,t2) .* data_var(t1,t2) ./ (data_var(t1,t2) + noise_sigma.^2);
    [t1,t2] = func_dyad2LH(j);
    ws(t1,t2) = wx(t1,t2) .* data_var(t1,t2) ./ (data_var(t1,t2) + noise_sigma.^2);
end

x_out  = func_IWT2_PO(ws, L, qmf);

%-------------------------------------------------------------------------
%------------------------------Inner Function ----------------------------
%-------------------------------------------------------------------------
% Estimate the signal variance value, which will be used for shrinkage
function result = func_signal_variance_estimation(wx, noise_sigma, win_size,ant_search_range_row_min,ant_search_range_row_max,ant_search_range_col_min,ant_search_range_col_max,J,L)

% System setup
ant_move_step_within_iteration = 5; % the numbe of iterations?
total_iteration_num = 5;
search_clique_mode = 8;

wx = abs(wx);
[nrow, ncol] = size(wx);
wx_1D = func_2D_LexicoOrder(wx);

% initialization
h = zeros(size(wx));
h_1D = func_2D_LexicoOrder(h);
p = 1./(wx.*(wx~=0)+1.*(wx==0)).*(wx~=0)+(wx==0);
p_1D = func_2D_LexicoOrder(p);

%paramete setting
alpha = 1;
beta = 2;
rho = 0.1;
w = 0.6;
A = 5000;
B = 10;

%use one ant for each pixel position
ant_total_num = nrow*ncol;
ant_current_row = zeros(ant_total_num, 1); % record the location of ant
ant_current_col = zeros(ant_total_num, 1); % record the location of ant
ant_current_val = zeros(ant_total_num, 1); % record the location of ant
rr = (meshgrid(1:nrow))'; cc = meshgrid(1:ncol);
ant_current_row = rr(:);
ant_current_col = cc(:);
ant_current_val = wx_1D((ant_current_row-1).*ncol+ant_current_col);

ant_search_range_row_min = ant_search_range_row_min(:);
ant_search_range_row_min = padarray(ant_search_range_row_min, [0 search_clique_mode-1],'replicate','post');
ant_search_range_row_max = ant_search_range_row_max(:);
ant_search_range_row_max = padarray(ant_search_range_row_max, [0 search_clique_mode-1],'replicate','post');
ant_search_range_col_min = ant_search_range_col_min(:);
ant_search_range_col_min = padarray(ant_search_range_col_min, [0 search_clique_mode-1],'replicate','post');
ant_search_range_col_max = ant_search_range_col_max(:);
ant_search_range_col_max = padarray(ant_search_range_col_max, [0 search_clique_mode-1],'replicate','post');

for nIteration = 1: total_iteration_num

        ant_current_path_val_mean = zeros(ant_total_num,1);
        ant_current_path_val_mean = ant_current_val;

        for nMoveStep = 1: ant_move_step_within_iteration-1

            if search_clique_mode == 4
                ant_search_range_row = [ant_current_row-1, ant_current_row, ant_current_row+1, ant_current_row];
                ant_search_range_col = [ant_current_col, ant_current_col+1, ant_current_col, ant_current_col-1];
            elseif search_clique_mode == 8
                ant_search_range_row = [ant_current_row-1, ant_current_row-1, ant_current_row-1, ant_current_row, ant_current_row,ant_current_row+1, ant_current_row+1, ant_current_row+1];
                ant_search_range_col = [ant_current_col-1, ant_current_col, ant_current_col+1, ant_current_col-1, ant_current_col+1, ant_current_col-1, ant_current_col, ant_current_col+1];
            end

            ant_current_row_extend = padarray(ant_current_row, [0 search_clique_mode-1],'replicate','post');
            ant_current_col_extend = padarray(ant_current_col, [0 search_clique_mode-1],'replicate','post');
            ant_search_range_val = zeros(ant_total_num,search_clique_mode);

            %replace the positions our of the image's range
            temp = (ant_search_range_row>=ant_search_range_row_min) & (ant_search_range_row<=ant_search_range_row_max) & (ant_search_range_col>=ant_search_range_col_min) & (ant_search_range_col<=ant_search_range_col_max);
            ant_search_range_row = temp.*ant_search_range_row + (~temp).*ant_current_row_extend;
            ant_search_range_col = temp.*ant_search_range_col + (~temp).*ant_current_col_extend;

            ant_search_range_transit_prob_h = zeros(size(ant_search_range_val));
            ant_search_range_transit_prob_p = zeros(size(ant_search_range_val));

            for ii=1:search_clique_mode
                ant_search_range_val(:,ii) = wx_1D((ant_search_range_row(:,ii)-1).*ncol+ant_search_range_col(:,ii));

                temp = ant_search_range_val(:,ii);
                temp = abs(temp - ant_current_path_val_mean);
                temp = 1./(temp.*(temp~=0)+1.*(temp==0)) .* (temp~=0)+(temp==0);

                ant_search_range_transit_prob_h(:,ii) = temp;
                ant_search_range_transit_prob_p(:,ii) = p_1D((ant_search_range_row(:,ii)-1).*ncol+ant_search_range_col(:,ii));

            end
            temp = (ant_search_range_transit_prob_h.^alpha) .* (ant_search_range_transit_prob_p.^beta);

            temp_sum = sum(temp,2);
            temp_sum = padarray(temp_sum, [0 search_clique_mode-1],'replicate','post');

            ant_search_range_transit_prob = temp ./ temp_sum;

            % generate a random number to determine the next position.
            rand('state', sum(100*clock));
            temp = rand(ant_total_num,1);
            temp = padarray(temp, [0 search_clique_mode-1],'replicate','post');
            temp = cumsum(ant_search_range_transit_prob,2)>=temp;
            temp = padarray(temp, [0 1],'pre');
            temp = (diff(temp,1,2)==1);

            temp_row = (ant_search_range_row .* temp)';
            [ii, jj, vv] = find(temp_row);
            ant_next_row = vv;
            temp_col = (ant_search_range_col .* temp)';
            [ii, jj, vv] = find(temp_col);
            ant_next_col = vv;

            ant_current_row = ant_next_row;
            ant_current_col = ant_next_col;
            ant_current_val = wx_1D((ant_current_row-1).*ncol+ant_current_col);
            ant_current_path_val_mean = (ant_current_path_val_mean.*nMoveStep + ant_current_val)/(nMoveStep+1);

            %update p;
            rr = ant_current_row;
            cc = ant_current_col;
            p_1D((rr-1).*ncol+cc,1) = w*p_1D((rr-1).*ncol+cc,1) + 1./(A+B.*ant_current_path_val_mean);
            p = func_LexicoOrder_2D(p_1D, nrow, ncol);

        end % end of nMoveStep

        p = (1-rho).*p;
        p_1D = func_2D_LexicoOrder(p);

end % end of nIteration

clear h h_1D temp
clear ant_current_row ant_current_col ant_current_val
clear ant_current_path_val_mean
clear ant_search_range_row search_range_path_col ant_search_range_val
clear ant_search_range_transit_prob_h ant_search_range_transit_prob_v ant_search_range_transit_prob
clear ant_search_range_row_min ant_search_range_row_max ant_search_range_col_min ant_search_range_col_max

result = wx.^2;

for j=(J-1):-1:L
    [t1,t2] = func_dyad2HH(j);
    result(t1,t2) = func_determine_class_fcm(wx(t1,t2), p(t1,t2), win_size, noise_sigma);
    [t1,t2] = func_dyad2HL(j);
    result(t1,t2) = func_determine_class_fcm(wx(t1,t2), p(t1,t2), win_size, noise_sigma);
    [t1,t2] = func_dyad2LH(j);
    result(t1,t2) = func_determine_class_fcm(wx(t1,t2), p(t1,t2), win_size, noise_sigma);
end

% ******************************************************************************
% **************************Inner Function *************************************
% ******************************************************************************
% Bi-class classification algorithm using Fuzzy C-means
function result = func_determine_class_fcm(wx, p, win_size, noise_sigma)

p_1D = func_2D_LexicoOrder(p);
wx_1D = func_2D_LexicoOrder(wx);
[nrow, ncol] = size(wx);

nFeature = p_1D(:)./(max(max(p_1D)));
[center,U,obj_fcn] = fcm(nFeature, 2,[2.0 100 1e-5 0]);

idx = zeros(nrow:ncol,1);
if sum(sum(wx_1D(U(1,:) >= U(2,:)))) >= sum(sum(wx_1D(U(1,:) < U(2,:))))
    idx = U(1,:) >= U(2,:);
else
    idx = U(1,:) < U(2,:);
end
idx = func_LexicoOrder_2D(idx, nrow, ncol);

center_class = idx(:);
padnum = win_size;
A = padarray(wx, [padnum padnum], 'replicate', 'bot');
B = padarray(idx, [padnum padnum], 'replicate', 'bot');

A = im2col(A, [win_size*2+1 win_size*2+1],'sliding');
A = A.^2;
B = im2col(B, [win_size*2+1 win_size*2+1],'sliding');

C = padarray(center_class', [size(A,1)-1, 0], 'replicate', 'post');

result = sum(A,1) ./ (win_size*2+1)./(win_size*2+1);
result = reshape(result, [nrow, ncol]);
result = (result>=noise_sigma^2).*(result-noise_sigma^2);

%-------------------------------------------------------------------------
%------------------------------Inner Function ----------------------------
%-------------------------------------------------------------------------
% Convert data from 2D format to 1D format
function result = func_2D_LexicoOrder(x)
temp = x';
result = temp(:);

%-------------------------------------------------------------------------
%------------------------------Inner Function ----------------------------
%-------------------------------------------------------------------------
% Convert data from 1D format to 2D format
function result = func_LexicoOrder_2D(x, nRow, nColumn)
result = reshape(x, nColumn, nRow)';

%-------------------------------------------------------------------------
%------------------------------Inner Function ----------------------------
%-------------------------------------------------------------------------
function result=min2(f)
%calculate minimum of 2D matrix
result=min(min(f));

%-------------------------------------------------------------------------
%------------------------------Inner Function ----------------------------
%-------------------------------------------------------------------------
function result=max2(f)
%calculate maximum of 2D matrix
result=max(max(f));

%-------------------------------------------------------------------------
%------------------------------Inner Function ----------------------------
%-------------------------------------------------------------------------
% Calculate the PSNR performance to two images
function result = func_psnr_gray(f, g)

f = double(f);
g = double(g);
Q=255; MSE=0;
[M,N]=size(f);
h = f - g;
MSE = sum(sum(h.*h));
MSE=MSE/M/N;
result=10*log10(Q*Q/MSE);

author: Sable

repo: mcbench-benchmarks


text:
function [data,compress,gridlocs,node_sizes,voxel_sizes] = fns_read_recipdata(infile)
%
% This function reads the reciprocity data and the related
% information into the MATLAB workspace.

% $Copyright (C) 2010 by Hung Dang$

% Read in the reciprocity data
data = hdf5read(infile,'/recipdata/data');

% Read in other parameters
compress    = hdf5read(infile,'/recipdata/roicompress') + 1;
gridlocs    = hdf5read(infile,'/recipdata/gridlocs');
node_sizes  = hdf5read(infile,'/model/node_sizes');
voxel_sizes = hdf5read(infile,'/model/voxel_sizes');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [sys,x0,str,ts]=NL_PID_fhan_c(t,x,u,flag,r0,h0,r1,h1,h,c,r2,h2)
switch flag,
    case 0,
        [sys,x0,str,ts]=mdlInitializeSizes(h);
    case 2,
        sys=mdlUpdates(x,u,r0,r1,h0,h1,h);
    case 3,
        sys=mdlOutputs(x,c,r2,h2);
    case 4,
        sys=mdlGetTimeOfNextVarHit(t,h);
    case {1,9},
        sys=[];
    otherwise
        error(['Unhandled flag=',num2str(flag)]);
end
function [sys,x0,str,ts]=mdlInitializeSizes(h)
    sizes=simsizes;
    sizes.NumContStates=0;
    sizes.NumDiscStates=5;
    sizes.NumOutputs=2;
    sizes.NumInputs=2;
    sizes.DirFeedthrough=1;
    sizes.NumSampleTimes=1;
    sys=simsizes(sizes);
    x0=[0;0;0;0;0];
    str=[];
    ts=[h 0];
function sys=mdlUpdates(x,u,r0,r1,h0,h1,h)
    fh0=fhan(x(1)-u(1),x(2),r0,h0);
    x(1)=x(1)+h*x(2);
    x(2)=x(2)+h*fh0;
    fh1=fhan(x(3)-u(2),x(4),r1,h1);
    x(3)=x(3)+h*x(4);
    x(4)=x(4)+h*fh1;
    x(5)=x(5)+h*(x(1)-x(3));
    sys=x;
function sys=mdlOutputs(x,c,r2,h2)
    sys(1)=-fhan(x(1)-x(3),c*(x(2)-x(4)),r2,h2);
    sys(2)=x(1);
function sys=mdlGetTimeOfNextVarHit(t,h)
    sys=t+h;
function y=fhan(x1,x2,r,h)
d=r*h;
d0=h*d;
y=x1+h*x2;
a0=sqrt(d^2+8*r*abs(y));
if abs(y)>d0
    a=x2+(a0-d)*sign(y)/2;
else
    a=x2+y/h;
end

if abs(a)>d
    y=-r*sign(a);
else
    y=-r*a/d;
end

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
function title = p07_title ( )

%*****************************************************************************80
%
%% P07_TITLE returns the title of problem p07.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'f(x) = exp ( 2 x ) if x < 0.5, 0 otherwise';

  return
end


text:
function [components,cliques,CC] = k_clique(k,M)
% k-clique algorithm for detecting overlapping communities in a network
% as defined in the paper "Uncovering the overlapping
% community structure of complex networks in nature and society" -
% G. Palla, I. Derényi, I. Farkas, and T. Vicsek - Nature 435, 814–818 (2005)
%
% [X,Y,Z] = k_clique(k,A)
%
% Inputs:
% k - clique size
% A - adjacency matrix
%
% Outputs:
% X - detected communities
% Y - all cliques (i.e. complete subgraphs that are not parts of larger
% complete subgraphs)
% Z - k-clique matrix
%
% Author : Anh-Dung Nguyen
% Email : anh-dung.nguyen@isae.fr

% The adjacency matrix of the example network presented in the paper
% M = [1 1 0 0 0 0 0 0 0 1;
%     1 1 1 1 1 1 1 0 0 1;
%     0 1 1 1 0 0 1 0 0 0;
%     0 1 1 1 1 1 1 0 0 0;
%     0 1 0 1 1 1 1 1 0 0;
%     0 1 0 1 1 1 1 1 0 0;
%     0 1 1 1 1 1 1 1 1 1;
%     0 0 0 0 1 1 1 1 1 1;
%     0 0 0 0 0 0 1 1 1 1;
%     1 1 0 0 0 0 1 1 1 1];

nb_nodes = size(M,1); % number of nodes

% Find the largest possible clique size via the degree sequence:
% Let {d1,d2,...,dk} be the degree sequence of a graph. The largest
% possible clique size of the graph is the maximum value k such that
% dk >= k-1
degree_sequence = sort(sum(M,2) - 1,'descend');
max_s = 0;
for i = 1:length(degree_sequence)
    if degree_sequence(i) >= i - 1
        max_s = i;
    else
        break;
    end
end

cliques = cell(0);
% Find all s-size kliques in the graph
for s = max_s:-1:3
    M_aux = M;
    % Looping over nodes
    for n = 1:nb_nodes
        A = n; % Set of nodes all linked to each other
        B = setdiff(find(M_aux(n,:)==1),n); % Set of nodes that are linked to each node in A, but not necessarily to the nodes in B
        C = transfer_nodes(A,B,s,M_aux); % Enlarging A by transferring nodes from B
        if ~isempty(C)
            for i = size(C,1)
                cliques = [cliques;{C(i,:)}];
            end
        end
        M_aux(n,:) = 0; % Remove the processed node
        M_aux(:,n) = 0;
    end
end

% Generating the clique-clique overlap matrix
CC = zeros(length(cliques));
for c1 = 1:length(cliques)
    for c2 = c1:length(cliques)
        if c1==c2
            CC(c1,c2) = numel(cliques{c1});
        else
            CC(c1,c2) = numel(intersect(cliques{c1},cliques{c2}));
            CC(c2,c1) = CC(c1,c2);
        end
    end
end

% Extracting the k-clique matrix from the clique-clique overlap matrix
% Off-diagonal elements <= k-1 --> 0
% Diagonal elements <= k --> 0
CC(eye(size(CC))==1) = CC(eye(size(CC))==1) - k;
CC(eye(size(CC))~=1) = CC(eye(size(CC))~=1) - k + 1;
CC(CC >= 0) = 1;
CC(CC < 0) = 0;

% Extracting components (or k-clique communities) from the k-clique matrix
components = [];
for i = 1:length(cliques)
    linked_cliques = find(CC(i,:)==1);
    new_component = [];
    for j = 1:length(linked_cliques)
        new_component = union(new_component,cliques{linked_cliques(j)});
    end
    found = false;
    if ~isempty(new_component)
        for j = 1:length(components)
            if all(ismember(new_component,components{j}))
                found = true;
            end
        end
        if ~found
            components = [components; {new_component}];
        end
    end
end

    function R = transfer_nodes(S1,S2,clique_size,C)
        % Recursive function to transfer nodes from set B to set A (as
        % defined above)

        % Check if the union of S1 and S2 or S1 is inside an already found larger
        % clique
        found_s12 = false;
        found_s1 = false;
        for c = 1:length(cliques)
            for cc = 1:size(cliques{c},1)
                if all(ismember(S1,cliques{c}(cc,:)))
                    found_s1 = true;
                end
                if all(ismember(union(S1,S2),cliques{c}(cc,:)))
                    found_s12 = true;
                    break;
                end
            end
        end

        if found_s12 || (length(S1) ~= clique_size && isempty(S2))
            % If the union of the sets A and B can be included in an
            % already found (larger) clique, the recursion is stepped back
            % to check other possibilities
            R = [];
        elseif length(S1) == clique_size;
            % The size of A reaches s, a new clique is found
            if found_s1
                R = [];
            else
                R = S1;
            end
        else
            % Check the remaining possible combinations of the neighbors
            % indices
            if isempty(find(S2>=max(S1),1))
                R = [];
            else
                R = [];
                for w = find(S2>=max(S1),1):length(S2)
                    S2_aux = S2;
                    S1_aux = S1;
                    S1_aux = [S1_aux S2_aux(w)];
                    S2_aux = setdiff(S2_aux(C(S2(w),S2_aux)==1),S2_aux(w));
                    R = [R;transfer_nodes(S1_aux,S2_aux,clique_size,C)];
                end
            end
        end
    end
end

repo: cluster


text:
function [b rho ypred matrix xval] = find_spatial_hist_fast(x)
% Function to find the spatial histogram for IMAGE and then compute Corr
% and fit 3-parameter logistic. returns logistic fitting value.
[nrows ncolms] = size(x);
 im = x;
 nbins = 100;
 max_x = 25;
 ss_fact = 4;
 [hval bins] = hist(x(:),nbins);
 x = map_matrix_to_closest_vec(im,bins);
for j = 1:ss_fact:max_x
    matrix = zeros(nbins);

    for i = 1:nbins

        [r c] = find(x==bins(i));

        if(isempty(r)~=1)
           h  = zeros(1,nbins);

            ctemp = c-j;
            rtemp = r;rtemp(ctemp<1) = [];
            ctemp(ctemp<1) = [];
            ind = sub2ind(size(im),rtemp,ctemp);
            if(~isempty(ind))
                h = h+hist(im(ind),bins);
            end

            ctemp = c+j; rtemp = r;rtemp(ctemp>ncolms) = [];
            ctemp(ctemp>ncolms) = [];
            ind = sub2ind(size(im),rtemp,ctemp);
             if(~isempty(ind))
                h = h+hist(im(ind),bins);
            end

            rtemp = r-j;ctemp = c;ctemp(rtemp<1) = [];
            rtemp(rtemp<1) = [];
            ind = sub2ind(size(im),rtemp,ctemp);
            if(~isempty(ind))
                h = h+hist(im(ind),bins);
            end

            rtemp = r+j;
            ctemp = c;ctemp(rtemp>nrows) = [];
            rtemp(rtemp>nrows) = [];
            ind = sub2ind(size(im),rtemp,ctemp);
            if(~isempty(ind))
                h = h+hist(im(ind),bins);
            end

            rtemp = r+j;ctemp = c;ctemp(rtemp>nrows) = [];rtemp(rtemp>nrows) = [];
            ctemp = ctemp+j;rtemp(ctemp>ncolms) = []; ctemp(ctemp>ncolms) = [];
            ind = sub2ind(size(im),rtemp,ctemp);
            if(~isempty(ind))
                h = h+hist(im(ind),bins);
            end

            rtemp = r+j;ctemp = c;ctemp(rtemp>nrows) = [];rtemp(rtemp>nrows) = [];
            ctemp = ctemp-j; rtemp(ctemp<1) = [];ctemp(ctemp<1) = [];
            ind = sub2ind(size(im),rtemp,ctemp);
            if(~isempty(ind))
                h = h+hist(im(ind),bins);
            end

            rtemp = r-j;ctemp = c;ctemp(rtemp<1) = [];rtemp(rtemp<1) = [];
            ctemp = ctemp+j;rtemp(ctemp>ncolms) = []; ctemp(ctemp>ncolms) = [];
            ind = sub2ind(size(im),rtemp,ctemp);
             if(~isempty(ind))
                h = h+hist(im(ind),bins);
            end

            rtemp = r-j;ctemp = c; ctemp(rtemp<1) = []; rtemp(rtemp<1) = [];
            ctemp = ctemp-j;rtemp(ctemp<1) = []; ctemp(ctemp<1) = [];
            ind = sub2ind(size(im),rtemp,ctemp);
             if(~isempty(ind))
                h = h+hist(im(ind),bins);
            end

            matrix(i,:) = h;
        end

    end
    X = [bins]'; Y = X; matrix = matrix/sum(sum(matrix));
    px = sum(matrix,2); py = sum(matrix,1)';
    mu_x = sum(X.*px); mu_y = sum(Y.*py);
    sigma2_x = sum((X-mu_x).^2.*px);    sigma2_y = sum((Y-mu_x).^2.*py);
    [xx yy] = meshgrid(X,Y);
    rho(j) = (sum(sum(xx.*yy.*matrix))-mu_x.*mu_y)/(sqrt(sigma2_x)*sqrt(sigma2_y));

end

rho = [rho(1:ss_fact:end)];
xval = [1:ss_fact:max_x];
% b = nlinfit([0 1:4:100],rho,@fit_spatial_corr,ones(1,5));
pp = polyfit([1:ss_fact:max_x],rho,3);
%  ypred = fit_spatial_corr(b,1:4:100);
ypred = polyval(pp,[1:ss_fact:max_x]);
err = sum((ypred-rho).^2);
b = [pp err];
%
% figure
% plot(rho,'r','LineWidth',3); hold on
% plot(ypred,'k'); title(['Error:',num2str( sum((ypred-rho).^2))]);
% toc
% subplot(1,3,1)
% imagesc(log(matrix+1));colormap(gray);axis xy
% subplot(1,3,2)
% bar(bins,sum(matrix,1)); axis([1 255 0 max(sum(matrix,1))])
% subplot(1,3,3)
% bar(bins,sum(matrix,2));axis([1 255 0 max(sum(matrix,2))])

repo: sr-metric


text:
function [HDR] = sclose(HDR)
% SCLOSE closes the file with the handle HDR
% [HDR] = sclose(HDR)
%    HDR.FILE.status = -1 if file could not be closed.
%    HDR.FILE.status = 0 indicates the file has been closed.
%
% see also: SOPEN, SREAD, SSEEK, STELL, SCLOSE, SWRITE

% 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.

%	$Id: sclose.m 2566 2010-11-24 13:04:34Z schloegl $
%	(C) 1997-2005,2006,2008 by Alois Schloegl <a.schloegl@ieee.org>
%    	This is part of the BIOSIG-toolbox http://biosig.sf.net/

if (HDR.FILE.FID<0) || ~HDR.FILE.OPEN,
        HDR.FILE.status = -1;
        %fprintf(HDR.FILE.stderr,'Warning SCLOSE (%s): invalid handle\n',HDR.FileName);
end;

if HDR.FILE.OPEN >= 2,          % write-open of files
	% check file length - simple test for file integrity
        EndPos = ftell(HDR.FILE.FID);          % get file length
        % force file pointer to the end, otherwise Matlab 6.5 R13 on PCWIN
        status = fseek(HDR.FILE.FID,0,'eof');

        if strcmp(HDR.TYPE,'BKR');
                if HDR.NS<1,
                        fprintf(HDR.FILE.stderr,'Error SCLOSE BKR: number of channels (HDR.NS) must be larger than zero.\n');
                        return;
                end;
                if HDR.NRec<1,
                        fprintf(HDR.FILE.stderr,'Error SCLOSE BKR: number of blocks (HDR.NRec) must be larger than zero.\n');
                        return;
                end;
                % check file length and write Headerinfo.
                HDR.SPR = (EndPos-HDR.HeadLen)/(HDR.NRec*HDR.NS*2);
                if isnan(HDR.SPR), HDR.SPR=0; end;
                if HDR.FILE.OPEN==3;
                        if any(isnan([HDR.NRec,HDR.NS,HDR.SPR,HDR.DigMax,HDR.PhysMax,HDR.SampleRate])), 	% if any unknown, ...
                                fprintf(HDR.FILE.stderr,'Error SCLOSE BKR: some important header information is still undefined (i.e. NaN).\n');
                                fprintf(HDR.FILE.stderr,'\t HDR.NRec,HDR.NS,HDR.SPR,HDR.DigMax,HDR.PhysMax,HDR.SampleRate must be defined.\n');
                                fprintf(HDR.FILE.stderr,'\t Try again.\n');
                        end;

                        fclose(HDR.FILE.FID);
                        HDR.FILE.FID = fopen(HDR.FileName,'r+');

                        count=fwrite(HDR.FILE.FID,HDR.VERSION,'int16');	        % version number of header
                        count=fwrite(HDR.FILE.FID,HDR.NS,'int16');	        % number of channels
                        count=fwrite(HDR.FILE.FID,HDR.SampleRate,'int16');      % sampling rate
                        count=fwrite(HDR.FILE.FID,HDR.NRec,'int32');            % number of trials: 1 for untriggered data
                        count=fwrite(HDR.FILE.FID,HDR.SPR,'uint32');            % samples/trial/channel
                        count=fwrite(HDR.FILE.FID,HDR.PhysMax,'int16');		% Kalibrierspannung
                        count=fwrite(HDR.FILE.FID,HDR.DigMax, 'int16');		% Kalibrierwert
                        count=fwrite(HDR.FILE.FID,zeros(4,1),'uint8');
                        count=fwrite(HDR.FILE.FID,[HDR.Filter.LowPass,HDR.Filter.HighPass],'float');

                        fseek(HDR.FILE.FID,32,'bof');
                        HDR.FLAG.TRIGGERED   = HDR.NRec>1;	% Trigger Flag
                        count = fwrite(HDR.FILE.FID,HDR.FLAG.TRIGGERED,'int16');           % FLAG TRIGGERED
                end;

	elseif strcmp(HDR.TYPE,'EDF') || strcmp(HDR.TYPE,'BDF') || strcmp(HDR.TYPE,'GDF'),
         	tmp = floor((EndPos - HDR.HeadLen) / HDR.AS.bpb);  % calculate number of records
                if isnan(tmp)
                	tmp = 0;
                end;

                if (HDR.NRec~=tmp)
                       if ~any(HDR.FILE.PERMISSION=='z')
                               HDR.NRec=tmp;
                               fseek(HDR.FILE.FID,236,'bof');
                               if strcmp(HDR.TYPE,'GDF')
                                        c=fwrite(HDR.FILE.FID,[HDR.NRec,0],'int32');
                               else
                                        fprintf(HDR.FILE.FID,'%-8i',HDR.NRec);
                                end;
                        else    %% due to a limitation in zlib
                                fprintf(HDR.FILE.stderr,'ERROR SCLOSE: number-of-records-field (HDR.NRec) could not be updated in file %s.\n',HDR.FileName);
                        end;
                end;

                if strcmp(HDR.TYPE,'GDF') && isfield(HDR,'EVENT'),
                	if ~all([HDR.NS, HDR.NRec, HDR.AS.bpb]>0)
                		HDR.AS.EVENTTABLEPOS = HDR.HeadLen;
                	else
	                        HDR.AS.EVENTTABLEPOS = HDR.HeadLen+HDR.AS.bpb*HDR.NRec;
	                end;
                        if isfield(HDR.EVENT,'VAL')
                                % handling of sparse (non-equidistant) sampling values
                                if isempty(HDR.EVENT.DUR)
                                        HDR.EVENT.DUR = zeros(size(HDR.EVENT.TYP));
                                end;
                                if (HDR.VERSION<1.90)
                                        fprintf(2,'Warning SCLOSE: GDF v2.0 or higher required for storing sparse samples but version is only %4.2f.\n',HDR.VERSION);
                                        %ix0 = ~isnan(HDR.EVENT.VAL);
                                        %HDR.EVENT.POS = HDR.EVENT.POS(~ix0);
                                        %HDR.EVENT.TYP = HDR.EVENT.TYP(~ix0);
                                        %HDR.EVENT.CHN = HDR.EVENT.CHN(~ix0);
                                        %HDR.EVENT.DUR = HDR.EVENT.DUR(~ix0);
                                        %HDR.EVENT = rmfield(HDR.EVENT,'VAL');
                                end;
	                end;
                        if isfield(HDR.EVENT,'VAL')
                                ix  = find(~isnan(HDR.EVENT.VAL));
                                tmp = unique(HDR.EVENT.CHN(ix));

                                flag.invalid = 0;
                                if (~all(HDR.EVENT.CHN(ix) > 0))
                                        flag.invalid = 1;
                                        fprintf(2,'Warning SCLOSE: Sparse sampling values without valid channel. Samples are not stored.\n');
                                        tmp = tmp(tmp > 0);
                                end;
                                if any(HDR.GDFTYP(tmp) > 6)
                                        fprintf(2,'Warning SCLOSE: Sparse sampling values must be of type integer using no more than 32 bits. Data is converted to uint32.\n');
                                end;

                                if any(HDR.AS.SPR(tmp)) || flag.invalid,
                                        fprintf(2,'Warning SCLOSE: Sparse sampling value for non-sparse channels not allowed. The following channel(s) is(are) affected: ');
                                        fprintf(2,'%i', tmp(HDR.AS.SPR(tmp)>0) );
                                        fprintf(2,'.  Samples are not stored.\n')
                                        ix0 = ~isnan(HDR.EVENT.VAL) & ~(HDR.EVENT.CHN>0); % initialize index with logical(0)
                                        for k = 1:length(tmp),
                                                if (HDR.AS.SPR(tmp(k))>0),
                                                        ix0 = ix0 | (HDR.EVENT.CHN==tmp(k));
                                                end;
                                        end;
                                        HDR.EVENT.POS = HDR.EVENT.POS(~ix0);
                                        HDR.EVENT.TYP = HDR.EVENT.TYP(~ix0);
                                        HDR.EVENT.CHN = HDR.EVENT.CHN(~ix0);
                                        HDR.EVENT.DUR = HDR.EVENT.DUR(~ix0);
                                        HDR.EVENT.VAL = HDR.EVENT.VAL(~ix0);
                                        ix = ~isnan(HDR.EVENT.VAL);
                                end;

                                % prepare for storing
                                HDR.EVENT.TYP(ix) = hex2dec('7fff');
                                HDR.EVENT.DUR(ix) = HDR.EVENT.VAL(ix);
                        end;

                        len = [length(HDR.EVENT.POS),length(HDR.EVENT.TYP)];
                        EVENT.Version = 1;
                        if isfield(HDR.EVENT,'CHN') && isfield(HDR.EVENT,'DUR'),
                                if any(HDR.EVENT.CHN) || any(HDR.EVENT.DUR),
                                        EVENT.Version = 3;
					len = [len,length(HDR.EVENT.CHN),length(HDR.EVENT.DUR)];
                                end;
                        end;

                        if any(len~=len(1))
                                fprintf(HDR.FILE.stderr,'Error SCLOSE-GDF: cannot write Event table, file %s not closed.\n',HDR.FileName);
				return;
                        else
                                status  = fseek(HDR.FILE.FID,HDR.HeadLen+HDR.AS.bpb*HDR.NRec,'bof');
                                %status = fseek(HDR.FILE.FID,0,'eof');
                                if ftell(HDR.FILE.FID)~=HDR.AS.EVENTTABLEPOS,
                                        fprintf(HDR.FILE.stderr,'Warning SCLOSE-GDF: inconsistent GDF file\n');
                                end
				if ~isfield(HDR.EVENT,'SampleRate'),
					HDR.EVENT.SampleRate = HDR.SampleRate;
				end;

				% write eventtable info: Version, SampleRate & length
				if (HDR.VERSION<1.94)
					tmp = HDR.EVENT.SampleRate;
					tmp = [EVENT.Version,mod(tmp,256),floor(mod(tmp,65536)/256),floor(tmp/65536)];
	                                fwrite(HDR.FILE.FID,tmp,'uint8');
	                                fwrite(HDR.FILE.FID,length(HDR.EVENT.POS),'uint32');
				else
					tmp = length(HDR.EVENT.POS);
					tmp = [EVENT.Version,mod(tmp,256),floor(mod(tmp,65536)/256),floor(tmp/65536)];
	                                fwrite(HDR.FILE.FID,tmp,'uint8');
	                                fwrite(HDR.FILE.FID,HDR.EVENT.SampleRate,'float32');
				end;

				% write table
                                c1 = fwrite(HDR.FILE.FID,HDR.EVENT.POS,'uint32');
                                c2 = fwrite(HDR.FILE.FID,HDR.EVENT.TYP,'uint16');
                                c3 = length(HDR.EVENT.POS); c4 = c3;
                                if EVENT.Version==3;
                                        c3 = fwrite(HDR.FILE.FID,HDR.EVENT.CHN,'uint16');
                                        c4 = fwrite(HDR.FILE.FID,HDR.EVENT.DUR,'uint32');
                                end;
                                if any([c1,c2,c3,c4]~=length(HDR.EVENT.POS))
                                        fprintf(2,'\nError SCLOSE: writing of EVENTTABLE failed. File %s not closed.\n', HDR.FileName);
                                        return;
                                end
                        end;
                end;

        elseif strcmp(HDR.TYPE,'CFWB');
                tmp = (EndPos-HDR.HeadLen)/HDR.AS.bpb;
                if isnan(tmp), tmp=0; end;
                if (tmp~=HDR.NRec);
                        if ~any(HDR.FILE.PERMISSION=='z')
                                HDR.NRec = tmp;
                                fseek(HDR.FILE.FID,14*4,-1);
                                count = fwrite(HDR.FILE.FID,HDR.NRec,'int32');           % channels
                        else
                                fprintf(HDR.FILE.stderr,'ERROR SCLOSE (CFWB): number-of-samples-field (HDR.SPR) could not be updated in file %s.\n',HDR.FileName);
                        end;
                end;

        elseif strcmp(HDR.TYPE,'SND');
                tmp = (EndPos-HDR.HeadLen)/HDR.AS.bpb;
                if isnan(tmp), tmp=0; end;
                if (HDR.FILE.OPEN==3) && (tmp~=HDR.SPR);
                        if ~any(HDR.FILE.PERMISSION=='z')
                                HDR.SPR = tmp;
                                fseek(HDR.FILE.FID,8,'bof');
                                count = fwrite(HDR.FILE.FID,HDR.SPR*HDR.AS.bpb,'uint32');           % bytes
                                fseek(HDR.FILE.FID,20,-1);
                                count = fwrite(HDR.FILE.FID,HDR.NS,'uint32');           % channels
                        else
                                fprintf(HDR.FILE.stderr,'ERROR SCLOSE (SND): fields HDR.SPR and HDR.NS could not be updated in file %s.\n',HDR.FileName);
                        end;
		end;

        elseif strcmp(HDR.TYPE,'AIF');
                if HDR.FILE.OPEN==3;
                        fseek(HDR.FILE.FID,4,-1);
                        count = fwrite(HDR.FILE.FID,EndPos-8,'uint32');           % bytes
                        fseek(HDR.FILE.FID,HDR.WAV.posis(2),-1);
                        count = fwrite(HDR.FILE.FID,EndPos-4-HDR.WAV.posis(2),'uint32');           % channels
                end;

        elseif strcmp(HDR.TYPE,'FLT') ;
		% warning FLT/header might have changed in needs to be rewritten
        elseif strcmp(HDR.TYPE,'WAV') ;
                if HDR.FILE.OPEN==3;
                        fseek(HDR.FILE.FID,4,-1);
                        count = fwrite(HDR.FILE.FID,EndPos-16,'uint32');           % bytes
                        fseek(HDR.FILE.FID,HDR.WAV.posis(2),-1);
                        count = fwrite(HDR.FILE.FID,EndPos-4-HDR.WAV.posis(2),'uint32');           % channels
                end;
        end;
end;

if 0,
elseif strcmp(HDR.TYPE,'ZIP')
        [SUCCESS,MESSAGE,MESSAGEID] = rmdir(HDR.ZIP.TEMPDIR,'s');

elseif strcmp(HDR.TYPE,'FIF') && HDR.FILE.OPEN;
        global FLAG_NUMBER_OF_OPEN_FIF_FILES
        rawdata('close');
        HDR.FILE.OPEN = 0;
        HDR.FILE.status = 0;
        FLAG_NUMBER_OF_OPEN_FIF_FILES = FLAG_NUMBER_OF_OPEN_FIF_FILES-1;

elseif strmatch(HDR.TYPE,{'ADI','GTEC','LABVIEW','MAT','MAT4','MAT5','native','XML','XML-FDA','SierraECG'}),
        HDR.FILE.OPEN = 0;

elseif HDR.FILE.OPEN,
        HDR.FILE.OPEN = 0;
        HDR.FILE.status = fclose(HDR.FILE.FID);
end;

author: PatternRecognition

repo: OpenBMI


text:
function [err, bf] = spm_bf (TR)
%
%   [err,] = spm_bf ()
%
%Purpose:
%
%
%
%Input Parameters:
%
%
%
%Output Parameters:
%   err : 0 No Problem
%       : 1  Problems
%
%
%
%Key Terms:
%
%More Info :
%
%
%
%
%     Author : Gang Chen
%     Date : Mon Oct 27 17:13:26 EST 2003
%     SSCC/NIMH/ National Institutes of Health, Bethesda MD 20892

%Define the function name for easy referencing
FuncName = 'spm_bf.m';

%Debug Flag
DBG = 1;

%initailize return variables
err = 1;

%	p(1) - delay of response (relative to onset)	   6        -gamd  in waver -GAM
%	p(2) - delay of undershoot (relative to onset)    16
%	p(3) - dispersion of response			   1
%	p(4) - dispersion of undershoot			   1   X
%	p(5) - ratio of response to undershoot		   6   X
%	p(6) - onset (seconds)				   0     -gamd
%	p(7) - length of kernel (seconds)		  32

      p      = [6 16 1 1 6 0 32];   %default parameters
	   [bf1 p]         = spm_hrf(TR,p);      %get gamma hrf from SPM function spm_hrf
	   dp     = 1;
	   p(6)   = p(6) + dp;    %for calculating time derivative
	   bf2    = (bf1 - spm_hrf(TR,p))/dp;
%		D      = (bf(:,1) - spm_hrf(TR,p))/dp;   %time derivative
		bf     = [bf1 bf2];     %combine the two

err = 0;
return;


text:
function r = subSet(r,ind)
% indexing of rotation
%
% Syntax
%   subSet(q,ind) %
%

r = subSet@quaternion(r,ind);
r.i = r.i(ind);

author: mtex-toolbox

repo: mtex

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


text:
function [lab]=vadlabel(data)

l=abs(data)>0;
n=find(l==0);
for i=2:length(n)-1
if data(n(i)-1)~=0 || data(n(i)+1)~=0
l(n(i))=1;
end
end
lab = l;


text:
clear all, close all, clc
addpath('./utils');
load allFaces.mat
X = faces(:,1:nfaces(1));
[L,S] = RPCA(X);

%%
inds = [3 4 14 15 17 18 19 20 21 32 43];
for k=[3 4 14 15 17 18 19 20 21 32 43]
    k
    subplot(2,2,1)
    imagesc(reshape(X(:,k),192,168)), colormap gray
    subplot(2,2,3)
    imagesc(reshape(L(:,k),192,168)), colormap gray
    subplot(2,2,4)
    imagesc(reshape(S(:,k),192,168)), colormap gray
    pause
end

author: dynamicslab


text:
function M = modelDiagnostics2(finput,varargin)
% M = modelDiagnostics2(finput,optional arguments)
%
% M = modelDiagnostics(stimList,'stimList','all','ISI',1.2,'TR',2,'contrasts',[1 1 -1 -1;1 -1 1 -1],'screen')
%
%M =
%        stimlist: your original list of conditions
%           delta: delta matrix sampled at .01 seconds
%           model: model sampled at .01 seconds
%       modelatTR: model sampled at TR
%        smoothed: temporally smoothed model at TR
%        conmodel: predictors for contrasts
%       consmooth: temporally smoothed predictors for contrasts
%           power: power (abs(fft)) matrix for smoothed model
%        conPower: power matrix for smoothed contrasts
%          energy: sum of squared deviations from mean for each predictor
%       conEnergy: sum of squared deviations from mean for each contrast
%             eff: efficiency of design.  1/trace(xtxi)
%  resampleeffinf: efficiency loss due to resampling compared to TR = 1 s.  eff(TR=1) / eff(TR)
%    smootheffinf: loss of efficiency due to smoothing; eff(nosmooth) / eff(smooth)
%         varbeta: variance of each individual beta
%          coneff: efficiency of contrast set
%      convarbeta: efficiency of each contrast beta
%             vif: variance inflation factors for predictors
%          conVif: variance inflation factors for contrasts
%            cond: condition number of model; larger is worse, > 30 is bad
%         conCond: condition number for contrast set
%           colin: correlation coefficients between individual predictors
%        conColin: correlation coefficients between contrasts
%       cbalfreqs: counterbalancing frequencies up to 3rd order
%            freq: frequencies with which each stimulus type occurs
%
% options: inptType is 'stimList', 'sampledModel','plot' (plot only)
%          outptType is 'all', 'short'
%
%
% More examples:
% M = modelDiagnostics(M.stimlist,'stimList','all','ISI',.85,'TR',2,'contrasts',[1 1 -1 -1],'screen')
% M2 = modelDiagnostics(M,'plot')
%
% Tor Wager, 11/17/01

%------------------------------------------------------------------------------
% Set up arguments and default values
%------------------------------------------------------------------------------
print2screen = 1;
outptType = 'all';
isfilt = 'n';
doverb = 0;
xc = [];

for i = 1:length(varargin)
    if isstr(varargin{i})
        switch varargin{i}
        case 'ISI', ISI = varargin{i+1};
        case 'TR', TR = varargin{i+1};
        case 'HRF', HRF = varargin{i+1};
        case 'HPlength', HPlength = varargin{i+1};
        case 'LPsmooth',LPsmooth = varargin{i+1};
        case 'contrasts',contrasts = varargin{i+1};
		case 'saturate',satval = varargin{i+1};
		case 'contrastweights',contrastweights = varargin{i+1};
		case 'maxOrder',maxOrder = varargin{i+1};
		case 'xc',xc = varargin{i+1};
		case 'screen',print2screen = 1;
		case 'noplot',print2screen = 0;
		case 'brief',outptType = 'brief';
		case 'plot',outptType = 'plot';
		case 'verbose',doverb = 1;
        end
    end
end

if doverb,disp('Computing model descriptives.'),end

% -----------------------------------------------------------------------------------
% * determine output type and just plot if only plot is specified
% -----------------------------------------------------------------------------------

switch outptType
case 'all'
    dopower = 1;doenergy = 1;doeff = 1;dovif = 1;docondnum = 1; dohrf = 1;,
    docolin = 1; docbal = 1;dofreq = 1;dotimes = 0;
case 'brief'
    if doverb,disp('Brief output selected.'),end
    dopower = 0; doenergy = 0; doeff = 0; dovif = 0; docondnum = 0; docolin = 1;
    docbal = 0; dofreq = 1; dotimes = 0; dohrf = 0;
case 'plot'
	M = finput;,if doverb,disp('Using already-computed M and plotting.'),end

	if ~isfield(M.ga,'TR'),M.ga.TR = input('M.ga.TR not found: Enter TR in s: ');,end
	if ~isfield(M.ga,'HPlength'),M.ga.HPlength = input('M.ga.HPlength not found: Enter HPlength in s: ');,end

	% plot predictors
    plotvarious(M,M.ga.TR,'Predictors')

	% make Q structure with contrast information
    Q.model = M.model * M.contrasts';
    try,Q.smoothed = M.consmooth;,catch,disp('Error: contrasts should be empty but they''re not?'),contrasts,end
	Q.modelatTR = resample(Q.model,1,M.ga.TR*10);

	if isempty(M.ga.HPlength), specstring = 'none';,else specstring = 'specify';,end
	[S,KL,KH] = use_spm_filter(M.ga.TR,size(Q.modelatTR,1),'hrf',specstring,M.ga.HPlength);
	if isempty(KH),KH = eye(size(Q.modelatTR));,end
    Q.filtered = Q.modelatTR - KH*(KH'*Q.modelatTR);
    Q.delta = M.delta * M.contrasts(:,1:end-1)';

	% plot contrasts
	plotvarious(Q,M.ga.TR,'Contrasts')

    return

otherwise eval(['help modelDiagnostics']),error('unknown output type')
end         % end switch

% 1. get every finput type to the stage of model at TR

% -----------------------------------------------------------------------------------
% * determine input type and set up.
% -----------------------------------------------------------------------------------
if 		isstruct(finput),inptType = 'Mstruct';
elseif 	size(finput,2) == 1, inptType = 'stimList';
elseif 	size(finput,2) > 1, inptType = 'model'
else	error('Illegal input.')
end

switch inptType
case 'Mstruct'
	if isfield(finput,'ga'),opt = finput.ga;
	else opt = finput;
	end

	M = finput;

    if isfield(M,'contrasts'), contrasts = M.contrasts;, end

    N = fieldnames(opt);
    for i = 1:length(N)
        eval([N{i} ' = opt.' N{i} ';'])
    end

	if ~isfield(M,'stimlist'), error('M must have stimlist field.'), end
	%if isfield(opt,'ISI'), ISI = opt.ISI;, end
	%if isfield(opt,'TR'), TR = opt.TR;, end
	%if isfield(opt,'nonlinthreshold'), satval = opt.nonlinthreshold;, end
	%if isfield(opt,'HPlength'), HPlength = opt.HPlength;, end
	%if isfield(opt,'LPsmooth'), LPsmooth = opt.LPsmooth;, end
	%if isfield(opt,'contrastweights'), contrastweights = opt.contrastweights;, end
	%if isfield(opt,'maxOrder'), maxOrder = opt.maxOrder;, end
	%if isfield(opt,'xc'), xc = opt.xc;, end
	%if isfield(opt,'contrasts'), contrasts = opt.contrasts;, end

    M.stimlist = double(M.stimlist);
	numsamps = ceil(size(M.stimlist,1) * ISI ./ TR);

case 'model'
	M.model = finput;
	M.modelatTR = M.model;
	if ~(sum(M.model(:,size(M.model,2))) == size(M.model,1))
		disp('No intercept detected in last column:  Adding intercept.')
		M.model(:,size(M.model,2)+1) = 1;
		M.modelatTR(:,size(M.modelatTR,2)+1) = 1;
        isfilt = input('Is the model filtered? (y/n) ','s');
        if strcmp(isfilt,'y')
            M.smoothed = M.modelatTR;
        end
	end
    numsamps = size(M.modelatTR,1);
    dohrf = 0;      % can''t do hrf without delta function input
    dofreq = 0;
    docbal = 0;

case 'stimList'

	numsamps = ceil(size(finput,1) * ISI ./ TR);
	M.stimlist = finput;
	M.stimlist = double(M.stimlist);

end 	% end switch

% -----------------------------------------------------------------------------------
% * build X model matrix, if necessary
% -----------------------------------------------------------------------------------
if ~strcmp(inptType,'model')
	if ~exist('ISI') == 1, ISI = input('Enter ISI in s: ');,end
	if ~exist('TR') == 1, TR = input('Enter TR in s: ');,end
	if ~exist('satval') == 1, satval = 2; ,end    % satval = input('Enter saturation threshold in s: ');,end
	if ~exist('HRF') == 1,
    	HRF = spm_hrf(.1);
    	HRF = HRF / max(HRF);
	end

	[M.modelatTR,M.model,M.delta] = rna2model(M.stimlist,ISI,HRF,TR,numsamps,satval,[]);
end

% 2. set up contrast and smoothing matrices

% -----------------------------------------------------------------------------------
% * set up filtering/smoothing and contrasts, do filtering
% -----------------------------------------------------------------------------------
if ~exist('TR') == 1, TR = input('Enter TR in s: ');,end
if ~exist('HPlength') == 1, HPlength = input('Enter HPlength in s: ');,end
if ~exist('LPsmooth') == 1, LPsmooth = input('Enter LPsmooth value (1,0): ');,end
if ~exist('contrasts') == 1, contrasts = input('Enter contrast matrix: ');,end
if ~exist('contrastweights') == 1, contrasts, contrastweights = input('Enter contrast weight vector: ');,end

if doverb,fprintf(1,'setting up filtering...'),end

% get smoothing matrix
[S,Vi,svi] = getSmoothing(HPlength,LPsmooth,TR,numsamps,xc);

% set up contrast weights
if isempty(contrastweights),contrastweights = ones(1,size(contrasts,1));,end

% set up contrasts for testing individual predictors (contrasts empty)
if isempty(contrasts),
	contrasts = eye(size(M.model,2) - 1);
	contrasts(:,end+1) = 0;
end

% make sure contrasts are the right length
if size(M.model,2) > size(contrasts,2)
		% warning('Model is larger than contrasts - extra intercept?  Adding to contrasts...')
		contrasts = [contrasts zeros(size(contrasts,1),size(M.model,2) - size(contrasts,2))];
end

if size(M.model,2) < size(contrasts,2)
    warning('Model is smaller than contrasts - extra zero in contrasts?  Truncating...')
	contrasts = contrasts(:,1:size(M.model,2));
end

if isempty(S), S = 1;,end
if ~(strcmp(inptType,'model')) | strcmp(isfilt,'n')
    M.smoothed = S * M.modelatTR;
end

% -----------------------------------------------------------------------------------
% * assign other output values to M structure
% -----------------------------------------------------------------------------------

M.contrasts = contrasts;
M.ga.contrasts = contrasts;
M.ga.TR = TR;
if exist('ISI') == 1, M.ga.ISI = ISI;, end
M.ga.HPlength = HPlength;
M.ga.LPsmooth = LPsmooth;
M.ga.contrastweights = contrastweights;
M.ga.xc = xc;

% -----------------------------------------------------------------------------------
% * make contrast-space predictors
% -----------------------------------------------------------------------------------
if isempty(contrasts)
    contrasts = eye(size(M.model,2));
end

try
	M.conmodel = M.modelatTR * contrasts';
catch
		whos contrasts
		M
		warning('Contrasts are different size from modelatTR.  This should not happen.  Adding intercept 0 to contrast and trying again...')
		contrasts = [contrasts zeros(size(contrasts,1),1)];
		M.contrasts = contrasts;
end

M.conmodel = M.modelatTR * contrasts';
M.consmooth = M.smoothed * contrasts';

% -----------------------------------------------------------------------------------
% * calculate individual field values
% -----------------------------------------------------------------------------------

if doenergy
    M.energy = var(M.smoothed) * size(M.smoothed,1);
    M.conEnergy = var(M.consmooth) * size(M.consmooth,1);
end

if doeff                                                        		% calculate efficiency
    if doverb,fprintf(1,'eff...'),end
	xtxitx = pinv(M.smoothed);                                       	% inv(X'S'SX)*(SX)'; pseudoinv of (S*X)
	[M.eff_fitness,M.eff] = calcEfficiency(M.ga.contrastweights,contrasts,xtxitx,svi);
    M.eff = M.eff';
	M.se_contrasts = 1 ./ M.eff;

    % ========= calc inefficiency of sampling and smoothing ===============
	myeff = 1/trace(inv(M.smoothed'*M.smoothed));                 % brief form of model efficiency, not considering autocorrelation
    noresampeff = 1/trace(inv(M.model'*M.model));                 % what would have happened at TR = .1?
    nosmootheff = 1/trace(inv(M.modelatTR'*M.modelatTR));
    M.resample_eff_loss = noresampeff / nosmootheff;
    M.smoothing_eff_loss = nosmootheff / myeff;

end

if dohrf
    if doverb,fprintf(1,'DX...'),end
    delta = [];
    for i = 1:max(M.stimlist(:,1))
        delta(:,i) = (M.stimlist == i);
    end
    %[M.DX] = tor_make_deconv_mtx2(M.delta,round(12 / M.ga.TR),10 * M.ga.TR);
    [M.DX] = tor_make_deconv_mtx2(delta,round(12 / M.ga.TR),M.ga.TR / M.ga.ISI);
    if ~isempty(S), M.DX = S * M.DX;,end
    xtxitx = pinv(M.DX);                                       		% inv(X'S'SX)*(SX)'; pseudoinv of (S*X)
    [M.hrf_eff_fitness,M.hrf_eff] = calcEfficiency([],[],xtxitx,svi);
    M.hrf_eff_avg = mean(M.hrf_eff(1:end-1));
end

if dovif
    if doverb,fprintf(1,'VIFs...'),end
	warning off
    M.vif = getvif(M.smoothed);
    if ~isempty(M.consmooth)
        M.conVif = getvif(M.consmooth);
    end
	warning on
end

if docondnum                                                    % calculate condition numbers
    M.cond = cond(M.smoothed);
    if ~isempty(M.consmooth)
        M.conCond = cond(M.consmooth);
    end
end
warning off
if docolin                                                      %   calculate predictor colinearity
    if doverb,fprintf(1,'colin...'),end
    % sqrt of average squared correlation
      colin = corrcoef(M.smoothed);
      M.colin = colin;
      if ~isempty(M.consmooth)
        colin = corrcoef(M.consmooth);
        M.conColin = colin;
      end
end
warning on
% calculate counterbalancing
if docbal
    if doverb,fprintf(1,'cbal...'),end
	if ~exist('maxOrder') == 1, maxOrder = input('Enter max counterbalancing order: ');,end

	conditions = 1:max(M.stimlist);
	freqConditions = ones(1,size(conditions,2)) / size(conditions,2);
    [M.cBal M.cbalfreqs] = getCounterBal(M.stimlist, maxOrder,conditions,freqConditions);
    % M.cbalfreqs = round(M.cbalfreqs * 100);
end

if dofreq
    for i = 1:size(conditions,2)
         M.freq(i) = sum(M.stimlist == conditions(i)) / size(M.stimlist,1);
    end
end

% calculate average time between stimulus repetitions in each trial type
if isfield(M,'stimlist')
    for i = 1:max(M.stimlist),M.timeBtwn(i) = mean(diff(find(M.stimlist==i))) .* M.ga.ISI;, end
end

if print2screen
    format short g
    format compact

    plotvarious(M,TR,'Predictors')
    if ~isempty(M.consmooth)
        Q.model = M.model * contrasts';
		Q.modelatTR  = M.modelatTR * contrasts';
        Q.smoothed = M.consmooth;
        Q.delta = M.delta * contrasts(:,1:end-1)';
        plotvarious(Q,TR,'Contrasts')
    end
end

if doverb
    disp(' ')
    disp('Correlated columns:')
    [x,y] = find(abs(corrcoef(xX.X)) > .5 & abs(corrcoef(xX.X)) < 1); [x y]
end
return

% subfunctions
% ================================================================== %
function plotvarious(M,TR,textlbl)

    if strcmp(textlbl,'Contrasts'), contit = 1;, else, contit = 0;, end

    % image the model
    % --------------------------------------

   figure; set(gcf,'Color','w')
    imagesc(M.modelatTR);colormap(copper)
    numsamp = size(M.smoothed,1);

    if contit
        title('Optimized Contrasts Heatmap, after filtering','FontSize',14)
    else
        title('Optimized Model Heatmap, after filtering','FontSize',14)
    end

    if isfield(M,'delta')
    	numregs = size(M.delta,2);
	% subtract one b/c last one is assumed to be intercept? sometimes...not always.

    else
	numregs = size(M.smoothed,2);	% subtract one b/c last one is assumed to be intercept? sometimes...not always.
    end

    nummodel = size(M.model,1);
    numsecs = numsamp * TR;
    deltascale = .1;
    x = 1:deltascale:numsecs+1;
    x = x(1:nummodel);
    x2 = 1:TR:numsecs+1;
    x2 = x2(1:numsamp);
    ymin = min(min(M.smoothed));
    ymax = max(max(M.smoothed));

    % Plot all columns in model after resampling and smoothing
    % --------------------------------------

    figure;set(gcf,'Color','w')
    for i = 1:numregs
  	  subplot(numregs,1,i)
      hold on
      plot(x,M.model(:,i),'k','LineWidth',1.5)
      plot(x2,M.modelatTR(:,i),'b','LineWidth',1.5)
      plot(x2,M.smoothed(:,i),'r','LineWidth',1.5)
      if isfield(M,'delta')
        for j = 1:size(M.delta,1)
            if M.delta(j,i) == 1
                plot([j*deltascale j*deltascale],[ymin 0],'k')
            end
        end
      end
      set(gca,'Ylim',[ymin ymax]);
      %axis([0 100 ymin ymax])
      set(gca,'XColor',[0 0 1])
      set(gca,'YColor',[0 0 1])
    end
    subplot(numregs,1,1);
    if contit
        title('Optimized Model Contrasts and Onsets','FontSize',14)
    else
        title('Optimized Model Regressors and Onsets','FontSize',14)
    end

    legend({'high-resolution' 'sampled at TR' 'filtered'})
    %title(textlbl,'FontSize',14)
    %title([textlbl ' Green: high-resolution   Blue: resampled at TR   Red: sampled and HP/LP filtered (smoothed)'])
    drawnow

    % Plot first column at various stages
    % --------------------------------------

    figure;set(gcf,'Color','w')
        subplot(3,1,1)
        ymin = min(min(M.model));
        ymax = max(max(M.model));
        plot(x,M.model(:,1),'k','LineWidth',2)
        title([textlbl '(#1): Original model - high samp rate'])
        set(gca,'Ylim',[ymin ymax]);
        if contit
            title('Optimized Model - first Contrast','FontSize',14)
        else
            title('Optimized Model - first Regressor','FontSize',14)
        end

        %axis([0 100 -1 2])
        subplot(3,1,2)
        plot(x2,M.modelatTR(:,1),'k','LineWidth',2)
        title('Resampled at TR')
        set(gca,'Ylim',[ymin ymax]);
        %axis([0 100 -1 2])
        subplot(3,1,3)
        plot(x2,M.smoothed(:,1),'k','LineWidth',2)
        title('Filtered')
        set(gca,'Ylim',[ymin ymax]);
        % axis([0 100 -1 2])

    drawnow

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function x = mysqueeze(x)
% Like squeeze(x), but row vectors become columns.

% Written by Thomas P Minka

s = size(x);
dim = find(s==1);
if ~isempty(dim)
  s(dim) = [];
  % [1 1] handles special cases
  x = reshape(x,[s 1 1]);
end

author: andrewssobral

repo: lrslibrary


text:
function varargout = minus(varargin)
%-   Subtraction of two SPHEREFUN objects.
%
%   F - G subtracts G from F, where F and G are SPHEREFUN objects or scalars.
%
% See also SPHEREFUN/PLUS, SPHEREFUN/UMINUS.

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

[varargout{1:nargout}] = minus@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun