text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the following exercises:
%   A) Cart and pendulum.
%              Simulate a cart and pendulum system. Then, try to control it
%              using a P controller.
%
% Help: Function prototype
% [y, t] = runge_kutta(f, y0, [t0 tfinal], timestep)
% where f is the function being integrated as dy/dt = f(t, y).
% y0 are the initial conditions
% t0: initial time
% tfinal: final time.
% h: time step for the calculations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cart_and_pendulum_controlled()
close all;

%uncomment to execute each of the exercises
exerciseA()%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simulate a cart and pendulum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function exerciseA()
global L
L = 0.8; % m length of the pendulum
t0 = 0;
tfinal = 10;
[t, x] = runge_kutta(@fcart_and_pendulum, [0 0 0.01 0.01]', [t0 tfinal], 0.01);
x = x(:, 1:length(t));

% animate the system
close all
figure
for i=1:length(t)
    %cart position
    posx_cart = x(1,i);
    posy_cart = 0;

    theta = x(3,i);
    %pendulum position
    posx_pend = L*sin(theta) + posx_cart;
    posy_pend = L*cos(theta);
    plot(posx_cart, posy_cart, 'rs', 'MarkerSize', 30), hold on
    plot(posx_pend, posy_pend, 'bo', 'MarkerSize', 12)
    xlim([-5 5])
    ylim([-0.9 0.9])
    line([posx_cart posx_pend], [posy_cart posy_pend])
    pause(0.05);
    hold off
end
% plot results
figure
plot(t, x(1,:), 'r'), hold
plot(t, x(2,:), 'g')
plot(t, x(3,:), 'b')
plot(t, x(4,:), 'k')
legend('Position x (m)', 'Speed xd (m/s)', '\theta (rad)', '\thetad (rad/s)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper function to solve a constant acceleration on x axis on a ramp.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xd = fcart_and_pendulum(t, x)
% We must return the solution of
% [dx1/dt; dx2/dt]
global L
M = 5; % kg mass of the cart
m = 5; % kg mass of the pendulum (centered at the tip)
g = 9.81; %m/s^ 2

% Computing a simple P controller
ref = 0; % theta reference
theta = x(3);
error = ref-theta;
k = 1000;
% the control action!
F = k*error;
% if F > 20
%     F=20;
% elseif F<-20
%     F=-20
% end

% compute xdd and thetadd to simulate the system
num = F + m*sin(x(3))*(g*cos(x(3)) - L*x(4)^2);
den = M + m*(1-cos(x(3))^2);
xdd = num/den;

num = F + g*(M+m)*tan(x(3))-m*L*x(4)^2*sin(x(3));
den = L*((M+m)/cos(x(3)) - m*cos(x(3)));
thetadd = num/den;

xd(1) = x(2);
xd(2) = xdd;
xd(3) = x(4);
xd(4) = thetadd;
xd = xd(:);

repo: ARTE


text:
function beta = linearCompressibility(S,x)
% computes the linear compressibility of an elasticity tensor
%
% Description
%
% $$\beta(x) = S_{ijkk} x_i x_j$$
%
% Input
%  S - elastic @complianceTensor
%  x - list of @vector3d
%
% Output
%  beta - linear compressibility in directions v
%

% return a function if required
if nargin == 1 || isempty(x)
  beta = S2FunHarmonicSym.quadrature(@(x) linearCompressibility(S,x),'bandwidth',2,S.CS);
  return
end

% compute tensor product
beta = EinsteinSum(S,[-1 -2 -3 -3],x,-1,x,-2);

author: mtex-toolbox

repo: mtex

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


text:
function varargout = ExportSnirfGUI(varargin)

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @ExportSnirfGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @ExportSnirfGUI_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

function ExportSnirfGUI_OpeningFcn(hObject, eventdata, handles, varargin)

if ~any([isa(varargin{1}, 'RunClass'), isa(varargin{1}, 'SubjClass'), isa(varargin{1}, 'SessClass'), isa(varargin{1}, 'GroupClass')])
    close(handles.exportSnirfGUI)
else
   handles.currElem = varargin{1};
end

handles.textElementName.String = handles.currElem.name;

[handles.dataKeys, handles.dataValues] = handles.currElem.procStream.output.GetAllData();

set(handles.listboxSelectData, 'String', handles.dataKeys);
set(handles.listboxSelectData, 'Max', length(handles.dataKeys));  % Enable multiple selection

handles.output = hObject;
guidata(hObject, handles);

% uiwait(handles.exportSnirfGUI);

% --- Outputs from this function are returned to the command line.
function varargout = ExportSnirfGUI_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;

function buttonSave_Callback(hObject, eventdata, handles)
selected = handles.listboxSelectData.Value;
data = handles.dataValues(selected);

% Prepare exported Snirf file
output = SnirfClass();
output.data = horzcat(data{:});
output.metaDataTags = MetaDataTagsClass();  % Generate default mdt
output.probe = handles.currElem.GetProbe();

name = handles.currElem.name;
name(name == '/') = '_';
name(name == '\') = '_';
name = erase(name, '.snirf');

if length(selected) < 1
    return
elseif length(selected) == 1
    suffix = handles.dataKeys{selected};
else
    suffix = 'exp';
end
[file, path, indx] = uiputfile(['homerOutput/', name, '_', suffix, '.snirf']);
if file ~= 0
    wb = waitbar(0.2, string(['Saving to ', file, ', please wait...']));
    output.Save([path, file]);
    waitbar(0.8, wb);
    close(wb);
    close(handles.exportSnirfGUI)
else
    close(handles.exportSnirfGUI)
end

function buttonCancel_Callback(hObject, eventdata, handles)
close(handles.exportSnirfGUI)

author: BUNPC


text:
%%*******************************************************************
%% schurmat_sblk: compute Schur complement matrix corresponding to
%%                SDP blocks.
%%
%% symm = 0, HKM
%%      = 1, NT
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

   function schur = schurmat_sblk(blk,At,par,schur,p,X,Y);

   global  nnzschur nzlistschur

   iter = par.iter;
   smallblkdim = par.smallblkdim;

   if isempty(smallblkdim); smallblkdim = 50; end
   if (nargin == 7); symm = 0; else; symm = 1; Y = X; end;
   m = length(schur);
   pblk = blk(p,:);
   if (iter == 1)
      nnzschur(size(blk,1),1) = m*m;
      nzlistschur = cell(size(blk,1),1);
   end
%%
   if (max(pblk{2}) > smallblkdim) | (length(pblk{2}) <= 10)
      %%
      %% compute schur for matrices that are very sparse.
      %%
      m1 = size(At{p,1},2);
      if issparse(schur); schur = full(schur); end;
      J = min(m1, max(find(par.nzlistA{p,1} < inf))-1);
      if (J > 0)
         if issparse(X{p}) & ~issparse(Y{p}); X{p} = full(X{p}); end
         if ~issparse(X{p}) & issparse(Y{p}); Y{p} = full(Y{p}); end
         if (iter <= 3)
            [nnzschur(p),nzlisttmp] = mexschur(pblk,At{p,1},par.nzlistA{p,1},...
            par.nzlistA{p,2},par.permA(p,:),Y{p},X{p},J,symm,schur);
            if (nnzschur(p) == mexnnz(nzlisttmp))
               nzlistschur{p} = nzlisttmp;
            else
               nzlistschur{p} = [];
            end
         else
            if isempty(nzlistschur{p})
               mexschur(pblk,At{p,1},par.nzlistA{p,1},...
               par.nzlistA{p,2},par.permA(p,:),Y{p},X{p},J,symm,schur);
            else
               mexschur(pblk,At{p,1},par.nzlistA{p,1},...
               par.nzlistA{p,2},par.permA(p,:),Y{p},X{p},J,symm,schur,nzlistschur{p});
            end
         end
      end
      %%
      %% compute schur for matrices that are not so sparse or dense.
      %%
      if (m1 < m) %% for low rank constraints
         ss = [0, cumsum(pblk{3})];
         len = sum(pblk{3});
         dd = At{p,3};
         DD = spconvert([dd(:,2:4); len,len,0]);
         XVD = X{p}*At{p,2}*DD;
         YVD = Y{p}*At{p,2}*DD;
      end
      L = max(find(par.nzlistAsum{p,1} < inf)) -1;
      if (J < L)
         len = par.nzlistAsum{p,1}(J+1); list = par.nzlistAsum{p,2}(1:len,:);
      end
      if (m1 > 0)
         for k = J+1:m
            if (k<=m1)
               isspAk = par.isspA(p,k);
               Ak = mexsmat(blk,At,isspAk,p,k);
               if (k <= L)
                  idx1 = par.nzlistAsum{p,1}(k)+1; idx2 = par.nzlistAsum{p,1}(k+1);
                  list = [list; par.nzlistAsum{p,2}(idx1:idx2,:)];
                  list = sortrows(list,[2 1]);
                  tmp = Prod3(pblk,X{p},Ak,Y{p},symm,list);
               else
                  tmp = Prod3(pblk,X{p},Ak,Y{p},symm);
               end
            else %%--- for low rank constraints
               idx = [ss(k-m1)+1 :ss(k-m1+1)];
               tmp = XVD(:,idx)* (Y{p}*At{p,2}(:,idx))';
            end
            if (~symm)
               tmp = 0.5*(mexsvec(pblk,tmp) + mexsvec(pblk,tmp'));
            else
               tmp = mexsvec(pblk,tmp);
            end
            permk = par.permA(p,k);
            idx  = par.permA(p,1:min(k,m1));
            tmp2 = schur(idx,permk) + mexinprod(blk,At,tmp,min(k,m1),p);
            schur(idx,permk) = tmp2;
            schur(permk,idx) = tmp2';
         end
      end
      if (m1 < m) %% for low rank constraints
         m2 = m - m1;
         XVtmp = XVD'*At{p,2};
         YVtmp = At{p,2}'*YVD;
         for k = 1:m2
            idx0 = [ss(k)+1 : ss(k+1)];
            tmp = XVtmp(:,idx0) .* YVtmp(:,idx0);
            tmp = tmp*ones(length(idx0),1);
            tmp3 = schur(m1+[1:m2],m1+k) + mexqops(pblk{3},tmp,ones(length(tmp),1),1);
            schur(m1+[1:m2],m1+k) = tmp3;
         end
      end
   else
      %%--- for SDP block where each sub-block is small dimensional
      if issparse(X{p}) & ~issparse(Y{p}); Y{p} = sparse(Y{p}); end
      if ~issparse(X{p}) & issparse(Y{p}); X{p} = sparse(X{p}); end
      tmp = mexskron(pblk,X{p},Y{p});
      schurtmp = At{p,1}'*tmp*At{p,1};
      %% schurtmp = 0.5*(schurtmp + schurtmp');
      if (norm(par.permA(p,:)-[1:m]) > 0)
         Perm = spconvert([(1:m)', par.permA(p,:)', ones(m,1)]);
         schur = schur + Perm'*schurtmp*Perm;
      else
         schur = schur + schurtmp;
      end
   end
%%*******************************************************************

author: yu-jiang

repo: radpbook

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


text:
%SOLVEEIG Solve eigenvalue problem.
%
%   [ U, L ] = SOLVEEIG( PROB, VARARGIN ) Solves the eigenvalue
%   problem described in the PROB finite element struct with
%   homogenous boundary conditions, and returns the eigenvectors U and
%   eigenvalues L. Accepts the following property/value pairs.
%
%       Input       Value/{Default}              Description
%       -----------------------------------------------------------------------------------
%       neigs       scalar/{6}                   Number of eigenvalues
%       sigma       scalar/string/{sm}           Range or type of eigenvalues
%                                                  see the eigs function for options
%       imass       scalar/{4}                   Mass matrix lumping
%                                                  1 - Full mass matrix
%                                                  2 - row sum lumping
%                                                  3 - diagonal lumping
%                                                  4 - HRZ diagonal lumping
%       init        u0|{expr}/{0}                Initial value solution or expression
%                                                for linearization point
%       solcomp     {all dvars}                  Dep. variables/subdomains to solve for
%       icub        scalar/{auto}                Cubature rule/order used in assembly
%                                                Default 1+max(shape function order)
%       isymm       scalar/{0}                   Symmetrize BCs if applicable
%       waitbar     scalar/{0}                   Show waitbar
%       fid         scalar/{1}                   File identifier for output ([]=no output)
%
%   See also EIGS, SOLVESTAT, SOLVETIME

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function [r,xy,v,zp,dfe] = weighted_reg_old2(y,w,varz)
% Calculate weighted average using weighted linear least squares
%
% :Model:
%
%   z_i = 1*zpop + noise
%
% :Inputs:
%
%   **Y:**
%        data matrix (nsub x T)
%
%   **w:**
%        weights
%
%   **varz:**
%        variance of data at each time point (nsub x T)
%
% :Outputs:
%
%   **r:**
%        weighted correlation coeff across columns of y
%
%   **xy:**
%        weighted covariance matrix
%
%   **v:**
%        weighted variance estimates for each column of y
%
%   **zp:**
%        weighted mean of each column of y
%

[m,n] = size(y);

% computation stuff
invxwx = inv(X'*W*X);
pxw = invxwx * X' * W;

W = diag(w);                    % Weight matrix

X = repmat(1,m,1);              % Design matrix - 1 column of all ones to calculate average

zp = px*y;        % weighted population mean

e = y - repmat(zp,m,1);         % residuals

dfe_v = zeros(n,1);

R = eye(m) - X*px;    % residual inducing matrix

% Calculate effective degrees of freedom

for i=1:n,

    V = diag(varz(:,i));

    dfe_v(i) = (trace(R*V)^2)/trace(R*V*R*V);       % Satherwaite approximation

end;

dfe = mean(dfe_v);               % Calculate average df over all time points.

MSE = (e'*W*e)/dfe;             % Mean square error

xy = inv(X'*W*X)*MSE;           % Covariance matrix for zp;

xy = 0.5*(xy+xy');              % Remove rounding error

v = diag(xy);                   % Variance for zp

r = xy./sqrt(v*v');             % Correlation matrix for zp

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [Results,u_star2,t]=Arm_1DOF_LinearStateTransition_Fun(StartAngle,Stable,I,l_c)
% this suggests an initial solution for the pendulium problem
%Stable 1 hanging pendulium, -1 inverted
%StartAngle given in degrees

%system parameters
m=1; %mass of plumb, kg
g=9.806; %gravity, m/sec^2
% I and l_c are now set by input

% simumation parameters
X_0=[StartAngle(1),0].'*pi/180; %intial value
X_F=[0,0].'*pi/180; %final value
t0=0;
tf=1;
T=linspace(t0,tf,1000);
T_2=[T, T(end)+(0.0001:0.01:0.1)];

% I*theta_ddot = -m*g*l_c*sin(theta)+u
% regular pendulium
% linearizing about theta=0 for now
A = [0 1; -Stable*m*g*l_c/I, 0];
B=[0; 1/I];

% W_c_fun=@(t) expm(A*t(n))*B*(B.')*expm(A.'*t(n));
% Written as subroutine below
temp11= quad(@(timein) W_c_fun(timein,A,B,1,1),t0,tf);
temp21= quad(@(timein) W_c_fun(timein,A,B,2,1),t0,tf);
temp12= quad(@(timein) W_c_fun(timein,A,B,1,2),t0,tf);
temp22= quad(@(timein) W_c_fun(timein,A,B,2,2),t0,tf);
W_c= [temp11, temp12; temp21, temp22];

u_pre=-B.';
u_post=(W_c\(expm(tf*A)*X_0-X_F));
u_star=zeros(size(T_2));
for alice=1:length(T)
u_star(alice)=u_pre*expm(A.'*(tf-T(alice)))*u_post;
end

Dynamics=@(t,x) [x(2); -Stable*(m*g*l_c/I)*sin(x(1))+...
    (sign(interp1(T_2,u_star,t,'pchip'))*min([abs(interp1(T_2,u_star,t,'pchip')),7]))/I];
[t,Results]=ode45(Dynamics,T_2,X_0);
% Apply saturation
u_star2=sign(interp1(T_2,u_star,t,'pchip')).*...
    min(abs(interp1(T_2,u_star,t,'pchip')),7);
end

function temp = W_c_fun(t,A,B,WhichRow,WhichCol)
temp=zeros([1,length(t)]);
for n=1:length(t)
    temp1=expm(A*t(n))*B*(B.')*expm(A.'*t(n));
    temp(n)=temp1(WhichRow,WhichCol);
end
end

author: Sable

repo: mcbench-benchmarks


text:
classdef MimCineTool < MimTool
    % MimCineTool. A tool for interactively moving through slices
    %
    %     MimCineTool is a tool class used with GemCinePanel to allow the user
    %     to cine through an image using mouse controls.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Cine'
        Cursor = 'arrow'
        ToolTip = 'Cine tool. Drag mouse to cine through slices.'
        Tag = 'Cine'
        ShortcutKey = 'n'
    end

    properties (Access = private)
        Callback
        ImageParameters
        StartCoords
        StartKPosition
    end

    methods
        function obj = MimCineTool(image_parameters, callback)
            obj.Callback = callback;
            obj.ImageParameters = image_parameters;
        end

        function MouseDragged(obj, screen_coords, last_coords)
            if ~isempty(obj.StartCoords)
                [min_coords, max_coords] = obj.Callback.GetImageLimits;
                coords_offset = screen_coords -  obj.StartCoords;

                y_range = max_coords(2) - min_coords(1);
                y_relative_movement = coords_offset(2)/y_range;
                direction = sign(y_relative_movement);
                y_relative_movement = abs(y_relative_movement);
                y_relative_movement = 100*y_relative_movement;
                y_relative_movement = ceil(y_relative_movement);

                k_position = obj.StartKPosition - direction*y_relative_movement;
                obj.ImageParameters.SliceNumber(obj.ImageParameters.Orientation) = k_position;
            end
        end

        function MouseDown(obj, screen_coords)
            obj.StartCoords = screen_coords;
            obj.StartKPosition = obj.ImageParameters.SliceNumber(obj.ImageParameters.Orientation);
        end

        function NewOrientation(obj)
            obj.StartCoords = [];
            obj.StartKPosition = [];
        end

        function ImageChanged(obj)
            obj.StartCoords = [];
            obj.StartKPosition = [];
        end

        function Enter(obj)
            obj.StartCoords = [];
            obj.StartKPosition = [];
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function this = apply_deformation_field(this, fileNameDeformationField, varargin)
% Applies a previoulsy estimated deformation field
%
%   Y = MrImage()
%   deformation_field = MrImage;
%   Y.apply_deformation_field(deformation_field)
%
% This is a method of class MrImage.
%
% IN
% MrImage object containing the deformation field
%
% OUT
%
% EXAMPLE
%   apply_deformation_field
%
%   See also MrImage

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2014-11-10
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.
%

this.save(this.get_filename('prefix', 'raw'));
matlabbatch = this.get_matlabbatch('apply_transformation_field', ...
    fileNameDeformationField, varargin{:});
save(fullfile(this.parameters.save.path, 'matlabbatch.mat'), ...
            'matlabbatch');
spm_jobman('run', matlabbatch);
% clean up: move/delete processed spm files, load new data into matrix

this.finish_processing_step('apply_transformation_field');

author: translationalneuromodeling

repo: tapas


text:
classdef (Abstract, Hidden = true) Base < adi.common.Attribute & ...
        adi.common.DebugAttribute & adi.common.RxTx & ...
        matlabshared.libiio.base & matlab.system.mixin.CustomIcon

    %adi.AD9361.Base Class
    %   This class contains shared parameters and methods between TX and RX
    %   classes
    properties (Nontunable)
        %SamplesPerFrame Samples Per Frame
        %   Number of samples per frame, specified as an even positive
        %   integer from 2 to 16,777,216. Using values less than 3660 can
        %   yield poor performance.
        SamplesPerFrame = 2^15;
    end

    properties (Nontunable, Logical)
        %EnableCustomFilter Enable Custom Filter
        %   Enable use of custom filter file to set SamplingRate,
        %   RFBandwidth, and FIR in datapaths
        EnableCustomFilter = false;
    end

    properties (Nontunable)
        %CustomFilterFileName Custom Filter File Name
        %   Path to custom filter file created from filter wizard
        CustomFilterFileName = '';
    end

    properties (Abstract)
        %CenterFrequency Center Frequency
        %   RF center frequency, specified in Hz as a scalar. The
        %   default is 2.4e9.  This property is tunable.
        CenterFrequency
        %SamplingRate Sampling Rate
        %   Baseband sampling rate in Hz, specified as a scalar
        %   from 65105 to 61.44e6 samples per second.
        SamplingRate
        %RFBandwidth RF Bandwidth
        %   RF Bandwidth of front-end analog filter in Hz, specified as a
        %   scalar from 200 kHz to 56 MHz.
        RFBandwidth
    end

    properties(Nontunable, Hidden)
        Timeout = Inf;
        kernelBuffersCount = 2;
        dataTypeStr = 'int16';
        phyDevName = 'ad9361-phy';
        iioDevPHY
    end

    properties (Hidden, Constant)
        ComplexData = true;
    end

    methods
        %% Constructor
        function obj = Base(varargin)
            coder.allowpcode('plain');
            obj = obj@matlabshared.libiio.base(varargin{:});
        end
        % Destructor
        function delete(obj)
            teardownLibad9361(obj);
            delete@adi.common.RxTx(obj);
        end
        % Check SamplesPerFrame
        function set.SamplesPerFrame(obj, value)
            validateattributes( value, { 'double','single' }, ...
                { 'real', 'positive','scalar', 'finite', 'nonnan', 'nonempty','integer','>',0,'<=',2^20}, ...
                '', 'SamplesPerFrame');
            obj.SamplesPerFrame = value;
        end
        % Check EnableCustomFilter
        function set.EnableCustomFilter(obj, value)
            validateattributes( value, { 'logical' }, ...
                { }, ...
                '', 'EnableCustomFilter');
            obj.EnableCustomFilter = value;
        end
        % Check CustomFilterFileName
        function set.CustomFilterFileName(obj, value)
            validateattributes( value, { 'char' }, ...
                { }, ...
                '', 'CustomFilterFileName');
            obj.CustomFilterFileName = value;
            if obj.EnableCustomFilter && obj.ConnectedToDevice %#ok<MCSUP>
                writeFilterFile(obj);
            end
        end
    end

    %% API Functions
    methods (Hidden, Access = protected)

        function icon = getIconImpl(obj)
            icon = sprintf(['AD9361 ',obj.Type]);
        end

        function setupLibad9361(obj)
            libName = 'libad9361';
            ad9361wrapperh = 'ad9361-wrapper.h';
            ad9361h = 'ad9361.h';
            fp = fileparts(which(ad9361h));
            loadlibraryArgs = {ad9361wrapperh,'includepath',fp,'addheader',ad9361h};
            if ~libisloaded(libName)
                msgID = 'MATLAB:loadlibrary:StructTypeExists';
                warnStruct = warning('off',msgID);
                [~, ~] = loadlibrary(libName, loadlibraryArgs{:});
                warning(warnStruct);
            end
            obj.iioDevPHY = calllib('libiio', 'iio_context_find_device',obj.iioCtx,'ad9361-phy');
        end

        function teardownLibad9361(~)
            libName = 'libad9361';
            if libisloaded(libName)
                unloadlibrary(libName);
            end
        end

        function writeFilterFile(obj)
            fir_data_file = obj.CustomFilterFileName;
            fir_data_str = fileread(fir_data_file);
            obj.setAttributeRAW('voltage0','filter_fir_en','0',false);
            obj.setAttributeRAW('voltage0','filter_fir_en','0',true);
            obj.setDeviceAttributeRAW('filter_fir_config',fir_data_str);
            obj.setAttributeRAW('voltage0','filter_fir_en','1',true);
            obj.setAttributeRAW('voltage0','filter_fir_en','1',false);
        end

    end

    %% External Dependency Methods
    methods (Hidden, Static)

        function tf = isSupportedContext(bldCfg)
            tf = matlabshared.libiio.ExternalDependency.isSupportedContext(bldCfg);
        end

        function updateBuildInfo(buildInfo, bldCfg)
            % Call the matlabshared.libiio.method first
            matlabshared.libiio.ExternalDependency.updateBuildInfo(buildInfo, bldCfg);
        end

        function bName = getDescriptiveName(~)
            bName = 'AD9361';
        end

    end
end

author: analogdevicesinc


text:
function out = ctranspose2x2(in)

% compute ctranspose of multiple 2x2 matrices, input is 2x2xN

out = conj(in);
out(1,2,:,:) = conj(in(2,1,:,:));
out(2,1,:,:) = conj(in(1,2,:,:));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function makespectrogramthumbnails(spectrogramFilename, spectrogramFraction)
debug.printfunctionstack('>');

% figure 1 should be a large spectrogram with traces, cropped nicely. Now remove labels and maximise panels.

ax=get(gcf, 'Children');

% Work out number of channels
pos1 = get(ax(1),'position'); % top trace
pos2 = get(ax(2),'position'); % top sgram
height_per_channel = pos1(4) + pos2(4);
numchannels = 0.95 / height_per_channel;
numpanels = numchannels * 2;

new_height_per_channel = 1 / numchannels;

% Remove all axes, tickmarks, labels, and axis boxes and title from view
for c=1:numpanels
    set(ax(c), 'Visible', 'off')
end

% Move panels
for channelNum = 1:numchannels
    [spectrogramPosition, tracePosition] = iceweb.calculatePanelPositions(numchannels, numchannels - channelNum + 1, spectrogramFraction, 0.0, 0.0, 1, 1);
    set(ax(channelNum*2 - 1), 'position', tracePosition);
    set(ax(channelNum*2), 'position', spectrogramPosition);
end

% we need a name for the labelless large spectrogram
[tmppath, tmpbase, tmpext] = fileparts(spectrogramFilename);
tmpfile = sprintf('%s/%s_labelless%s',tmppath,tmpbase,tmpext);

% print large labelless PNG
iceweb.saveImageFile(tmpfile, 72);

% load then delete temporary file
I = imread(tmpfile);
delete(tmpfile)

% Resize the image (aspect ratio 16:21 same as 576:756) and convert it to an indexed image with 256 colors
% (Note: we were originally creating 150x96, which is far off the aspect ratio of large spectrograms)
%[X,map] = rgb2ind(imresize(I, [126 96]), 256);
%thumbnailfile = sprintf('%s/smallest_%s%s',tmppath, tmpbase, tmpext);
%imwrite(X,map,thumbnailfile,'PNG');
%[X,map] = rgb2ind(imresize(I, [147 112]), 256);
%thumbnailfile = sprintf('%s/smaller_%s%s',tmppath, tmpbase, tmpext);
%imwrite(X,map,thumbnailfile,'PNG');
[X,map] = rgb2ind(imresize(I, [198 151]), 256);
thumbnailfile = sprintf('%s/%s_thumb%s',tmppath, tmpbase, tmpext);
imwrite(X,map,thumbnailfile,'PNG');
close;

debug.printfunctionstack('<');

author: geoscience-community-codes

repo: GISMO


text:
function f = make(varargin)
%MAKE   Constructor shortcut for TRIGTECH objects.
%   For TRIGTECH methods implemented at the TRIGTECH level, it is not possible
%   to call the class constructor file corresponding to a TRIGTECH object
%   directly. F = MAKE(VARARGIN) allows us to get around this and construct a
%   TRIGTECH.
%
% See also TRIGTECH.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Developer note: This is essentially a "factory method" in the sense of
% _Design Patterns_ by Gamma, Helm, Johnson, and Vlissides.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f = trigtech(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


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

disp('Installing the HDR Toolbox...');

folder = cellstr('Alignment');
folder = cellstr('Analysis');
folder = [folder, cellstr('Alignment/util')];
folder = [folder, cellstr('BatchFunctions')];
folder = [folder, cellstr('ColorCorrection')];
folder = [folder, cellstr('ColorSpace')];
folder = [folder, cellstr('Compression')];
folder = [folder, cellstr('Deghosting')];
folder = [folder, cellstr('Deghosting/util')];
folder = [folder, cellstr('EnvironmentMaps')];
folder = [folder, cellstr('EO')];
folder = [folder, cellstr('EO/util')];
folder = [folder, cellstr('Formats')];
folder = [folder, cellstr('Generation')];
folder = [folder, cellstr('Generation/util')];
folder = [folder, cellstr('Generation_video')];
folder = [folder, cellstr('Generation_video/util')];
folder = [folder, cellstr('IBL')];
folder = [folder, cellstr('IBL/util')];
folder = [folder, cellstr('IO')];
folder = [folder, cellstr('IO_stack')];
folder = [folder, cellstr('IO_video')];
folder = [folder, cellstr('LaplacianPyramids')];
folder = [folder, cellstr('Metrics')];
folder = [folder, cellstr('Metrics/util')];
folder = [folder, cellstr('NativeVisualization')];
folder = [folder, cellstr('Tmo')];
folder = [folder, cellstr('Tmo/util')];
folder = [folder, cellstr('Tmo_video')];
folder = [folder, cellstr('Tools')];
folder = [folder, cellstr('util')];

for i=1:length(folder)
    addpath([pwd(), '/source_code/', char(folder(i))], '-begin');
end

addpath([pwd(), '/demos/'], '-begin');
addpath([pwd(), '/demos/ldr_v'], '-begin');
addpath([pwd(), '/demos/output'], '-begin');
addpath([pwd(), '/demos/stack'], '-begin');
addpath([pwd(), '/demos/stack_alignment'], '-begin');
addpath([pwd(), '/demos/stack_ghost'], '-begin');

savepath

clear('folder');

disp('done!');

disp(' ');
disp('Check demos in the folder ''demos'' for learning how to use the HDR Toolbox!');
disp(' ');
disp('If you use the toolbox in your research, please reference the book in your papers:');
disp('@book{Banterle:2017,');
disp('      author = {Banterle, Francesco and Artusi, Alessandro and Debattista, Kurt and Chalmers, Alan},');
disp('      title = {Advanced High Dynamic Range Imaging (2nd Edition)},');
disp('      year = {2017},');
disp('      month={July},');
disp('      isbn = {9781498706940},');
disp('      publisher = {AK Peters (CRC Press)},');
disp('      address = {Natick, MA, USA},');
disp('      }');

disp(' ');
disp('NOTE ON TONE MAPPING:');
disp('=====================');
disp('The majority of TMOs return tone mapped images with linear values. This means that gamma encoding');
disp('needs to be applied to the output of these TMOs before visualization or before writing tone mapped images');
disp(' on the disk; otherwise these images may appear dark.');
disp('A few operators (e.g. Mertens et al.''s operator) return gamma encoded values,');
disp('so there is no need to apply gamma to them; in this case a message (e.g. a Warning) is displayed');
disp('after tone mapping alerting that there is no need of gamma encoding.');

disp(' ');
disp('NOTE ON EXPANSION OPERATORS (INVERSE/REVERSE TONE MAPPING):');
disp('=====================');
disp('The majority of EOs require to have as input LDR images in the range [0,1] that are LINEARIZED.');
disp('This means that the camera response function (CRF) or the gamma encoding has been removed.');
disp('This operation is MANDATORY in order to generate fair comparisons; please do use the gammaRemoval');
disp('parameter to remove gamma if you do not have the CRF of the input image. RAW files do not require ');
disp('this step because they are already linear.');

str_cur = pwd();
try
    cd([str_cur, '/source_code/IO/']);
    mex('read_exr.cpp');
    mex('write_exr.cpp');
catch err
    disp('read_exr.cpp and write_exr.cpp were not compiled');
    mex -setup cpp;
    try
        mex('read_exr.cpp');
        mex('write_exr.cpp');
    catch err
        disp('read_exr.cpp and write_exr.cpp were not compiled');
    end
end

try
    cd([str_cur, '/source_code/util/']);
    mex('bilateralFilterS.cpp');
catch err
    disp('bilateralFilterS.cpp was not compiled');
    mex -setup cpp;
    try
        mex('bilateralFilterS.cpp');
    catch err
        disp('bilateralFilterS.cpp was not compiled');
    end
end

cd(str_cur);

author: banterle


text:
function [fBetaBstMean, fBetaBstStd, mResult] = calc_beta_bootstrap(mCat,nBst,nSample)
% function [fBetaBstMean, fBetaBstStd, mResult] = calc_beta_bootstrap(mCat,nBst,nSample)
% -------------------------------------------------------------------------
%
% Calculate beta uncertainty for Michaels` Stress Tensor inversion approach boostrapping
% the Focal Mechanisms of the earthquake catalog mCat. This is a stability
% test of the method. No flipping of focal mechanisms included.
%
% Incoming
% mCat: EQ catalog or [dip direction dip rake]
% nBst: Number of boostraps
% nSample: Number of events from catalog to be used, default is all
%

if ~exist('nSample','var')
    nSample = length(mCat(:,1));
end

% Bootstrap event to be used
[nNy,nNx] = size(mCat);
if nNx > 3
    mCat = mCat(:,10:12);
end

% Bootstrap numbers
vNy = 1:1:nNy;
mBootstrap = bootrsp(vNy,nBst);

%Path
sZmapPath = './AddOneFiles/zmap/';

% Container for bootstrap beta values
mResult = [];

for nCnt = 1:1:nBst
    mCatBst = mCat(mBootstrap(1:nSample,nCnt),:);
    % Array of focal mechanisms: dip direction, dip, rake
    mFPS = mCatBst;

    % Do inversion using A. Michael code
    % Create file for inversion
    fid = fopen('data2','w');
    str = ['Inversion data'];str = str';
    fprintf(fid,'%s  \n',str');
    fprintf(fid,'%7.3f  %7.3f  %7.3f\n',mFPS');
    fclose(fid);
    % slick calculates the best solution for the stress tensor according to
    % Michael(1987): creates data2.oput
    sPath = pwd;
    unix([sZmapPath 'external/slick ' sPath '/data2 ']);

    % Get data from data2.oput
    sFilename = ['data2.oput'];
    % Calculate average angle between tangential traction predicted by best
    % stress tensor and the slip direction
    [fBeta, fStdBeta, fTauFit, fAvgTau, fStdTau] = import_slickoput(sFilename);

    % Delete existing from earlier runs data2.slboot
    sData2 = ['data2.slboot'];
    delete(sData2);

    % Stress tensor inversion
    unix([sZmapPath 'external/slfast ' sPath '/data2 ']);
    sGetFile = ['data2.slboot'];
    load(sGetFile);
    % Description of data2
    % Line 1: Variance S11 S12 S13 S22 S23 S33 => Variance and components of
    % stress tensor (S = sigma)
    % Line 2: Phi S1t S1p S2t S2p S3t S3p
    % Phi is relative size S2/S1, t=trend, p=plunge (other description)
    fVariance = data2(1,1);
%     fS11 = data2(1,2);
%     fS12 = data2(1,3);
%     fS13 = data2(1,4);
%     fS22 = data2(1,5);
%     fS23 = data2(1,6);
%     fS33 = data2(1,7);
    fPhi = data2(2,1);
    fS1Trend = data2(2,2);
    fS1Plunge = data2(2,3);
    fS2Trend = data2(2,4);
    fS2Plunge = data2(2,5);
    fS3Trend = data2(2,6);
    fS3Plunge = data2(2,7);
    % Container
    mResult = [mResult; fBeta fStdBeta fTauFit fAvgTau fStdTau fVariance fPhi fS1Trend fS1Plunge fS2Trend fS3Trend fS3Plunge];
end

% Bootstrap standar deviation
fBetaBstStd = calc_StdDev(mResult(:,1));
fBetaBstMean = mean(mResult(:,1));

author: CelsoReyes


text:
function params = er_defaultParams
%
%  params = er_defaultParams;
%
% Returns default parameters for running
% event-related analyses. Chose them
% to be nice for most KGS-lab block-design
% extps.
%
%
% ras 04/05.
% ras 11/06: added HRF sub-params (for standard H.R.F.s)
% dar 3/07    fixed typo in glmHRF_params, added comments on each parameter
params.eventAnalysis = 1; % back compatibility w/ mrLoadRet code

% detrend flag:
%--------------
% -1 linear detrend, 0 no detrend, 1 multiple boxcar smoothing,
% 2 quartic trend removal
params.detrend = 1;

% detrend frames:
% for detrend option 1, # of frames for smoothing
% kernal (roughly equivalent to cutoff of highpass
% filter):
params.detrendFrames = 20;

% inhomogeneity correction flag:
%-------------------------------
%     0 do nothing
%     1 divide by the mean, independently at each voxel
%     2 divide by null condition
%     3 divide by spatial gradient estimate
params.inhomoCorrect = 1;

% temporal normalization flag:
params.temporalNormalization = 0;

% seconds relative to trial onset to take for each trial
params.timeWindow = -8:24;

% period to look for peaks in t-tests, in seconds
params.peakPeriod = 4:14;

% period to use as baseline in t-tests, in seconds
params.bslPeriod = -8:0;

% threshold for significant activations
params.alpha = 0.05;

% # secs to shift onsets in parfiles, relative to time course
params.onsetDelta = 0;

% conditions to use for calculating signal-to-noise, HRF
params.snrConds = 1;

% flag for which hemodynamic impulse response
% function to use if applying a GLM:
% -------------------------------------------
% 0: deconovolve (selective averaging)
% 1: estimate HRF from mean response to all non-null conditions
% 2: Use Boynton et all 1998 gamma function
% 3: Use SPM difference-of-gammas
% 4: Use HRF from Dale and Buckner, 1997 (very similar to Boynton
%    gamma)
% OR, if flag is a char: name of a saved HRF function
%    (stored in subject/HRFs/, where subject is the subject's
%     anatomy directory, where the vAnatomy is stored)
params.glmHRF = 2;                    % boynton hrf
params.glmHRF_params = [3 1.08 2.05]; % good params for Boynton HIRF

% flag for whether or not to estimate temporally-correlated
% noise in data when applying a GLM, referred to as 'whitening':
% (see Dale and Burock, HBM, 2000):
params.glmWhiten = 0;

% # of events per block: this indicates whether each event
% is actually a block of evenly-spaced sub-events (trials).
params.eventsPerBlock = 1; % if block-design, # of events per block

params.ampType = 'betas';

% flag to zero baseline or not
params.normBsl = 1;

% params.colorOrder = []; % will get from tc_colorOrder, but
%                         % depends on # of conditions specified
%                         % in parfiles.

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% Peak of spectrogram frequency estimation
%
% Estimates the instantaneous frequency of the input signal by extracting
% the peaks of the spectrogram.
%
%
% Usage:
%
%     ife = sfpe (signal, window_length, time_res [, fft_length] )
%
% Parameters:
%
%    signal
%
%	  Input one dimensional signal to be analysed.	An analytic signal
%	  is required for this function, however, if signal is real, a
%	  default analytic transformer routine will be called from this
%	  function before computing spectrogram.
%
%    window_length
%
%	  The size of the kernel used for analysis. window_length must
%	  be odd. The kernel used will be defined from -(window_length+1)/2 to
%	  +(window_length+1)/2 in both time and lag dimensions.
%
%    time_res
%
%	  The number of time samples to skip between successive slices of
%	  the analysis.
%
%    fft_length
%
%         Zero-padding at the fft stage of the analysis may be specified by
%         giving an fft_length larger than normal. If fft_length is not
%	  specified then the next highest power of two above the
%         signal length will be used. If fft_length is not a power
%         of two, then the next highest power of two is used.
%
%
%
%  See Also: spec
% TFSAP 7.0
% Copyright Prof. B. Boashash
% Qatar University, Doha
% email: tfsap.research@gmail.com

author: HeLiangHIT


text:
function [VecTot, TOT, timei] = ierm(P, Mat, NumDesc)

%ierm returns a matrix containing the best models for all the paths of the
%Enhanced Replacement Method using an initial set of descriptors with a
%very high S, given by an inverse RM.
%TOT contains all the relative results of each step showing the evolution of the method.
%
%
%	   Input:
%             P             Property vector
%             Mat           Descriptors matrix with descriptors pool
%             NumDesc       Number of descriptors that the model will have
%
%     Returns:
%
%            VecTot           vector containing the best model for all the
%                               paths of the Replacement Method
%            TOT               contains all the relative results
%                               showing the evolution of the method.
%
% Andrew G. Mercader, Pablo R. Duchowicz
% INIFTA, La Plata, Argentina
% Created: 12 Nov 2007

timei=cputime;

if (nargin < 3)
   error('the function requires at least 5 input variables. Type ''help ierm''.');
end

[c_m, r_m] = size(Mat);

NTot=r_m;

VecI=1:NumDesc;
lindep = ld(P, VecI,Mat);

if lindep==100
    for j=1:1000000000000000000000000000000000000
     VecI = randint(1,NumDesc,[1, NTot]);
     lindep = ld(P, VecI,Mat);
        if lindep==0
            break;
        end
    end
end

[VecT] = rmt_inv(P, VecI, Mat);
VecJ=VecT(NumDesc,3:end);

[VecTot, TOT, time] = erm(P, VecJ, Mat);

timei=cputime-timei
end

%End


text:
function [x,P]=Kfilter(v,H,R,x_pre,P_pre)

nstate=size(x_pre,1);
K=P_pre*H'*(H*P_pre*H'+R)^-1;
x=x_pre+K*v;
P=(eye(nstate)-K*H)*P_pre;

return

repo: GINav


text:
function f = flops(fl)
% FLOPS         Get or set the global flopcount variable.
% FLOPS returns the current flopcount.
% FLOPS(F) sets flopcount to F.
%
% 0 flops: -x ' repmat
% 1 flop each: + - .*
% 2 flops each: < > == ~=
% For complex numbers, + is 2 flops, * is 6 flops.
% col_sum(x) takes (rows(x)-1)*cols(x) flops (use FLOPS_COL_SUM).
% row_sum(x) takes rows(x)*(cols(x)-1) flops (use FLOPS_ROW_SUM).
% Use FLOPS_DIV for ./
% Use FLOPS_RANDNORM for randn
% Use FLOPS_SQRT for sqrt
% Use FLOPS_POW for .^
% Use FLOPS_EXP for exp
% Use FLOPS_LOG for log, sin, and other special functions.
%
% See FLOPS_MUL, FLOPS_SOLVE, FLOPS_INV, FLOPS_CHOL, FLOPS_DET, ...

global flopcount;
if nargin == 1
  flopcount = fl;
  if nargout == 1
    f = fl;
  end
else
  f = flopcount;
end

author: Cloud-CV

repo: object-proposals


text:
function [Min_numerator,r]=ReGet_min_numerator(x1,x2,N,m,r,fs)
%******************************************************
% $ This function is usded to increase r value if the 'at least 5 match' is not achieved.
% $ Used for both COSEn and COFMEn.
%
% $ Variable declaration:
% x1 and x2 are the vector at dimentions of m and m+1
% N is time series length
% m is embedding dimension (usually m=1)
% r is changed threshold value
% fs sample rate
%
% $ Author: Chengyu Liu (bestlcy@sdu.edu.cn)
%           Institute of Biomedical Engineering,
%           Shandong University
% $Last updated:  2015.10.10
% Las updated: 2017.19.12 (by Giulia Da Poian) vectorized
%
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information

dim = N-m;

% Fast implementation

t1 = repmat(x2(:,1),[1 dim]);
t2 = toeplitz([x2(1,1); flipud(x2(2:end,1))], x2(:,1));
t3 = repmat(x2(:,2),[1 dim]);
t4 = toeplitz([x2(1,2); flipud(x2(2:end,2))], x2(:,2));
d2max = max(abs(t1-t2),abs(t3-t4));
d2max = cell2mat(arrayfun(@(x) circshift(d2max(x,:),[1 dim-x+1]),(1:dim)','un',0));
Min_numerator = numel(find(d2max(dim,:)<=r));

% Old implementation (slower)
% % for i=1:N-m
% %     Min_numerator=0;
% %     for j=1:N-m
% %         d2max(i,j)=max(abs(x2(i,:)-x2(j,:)));
% %         if  d2max(i,j)<=r
% %             Min_numerator=Min_numerator+1
% %         end
% %     end
% % end

if N<20
    r=r+1;
else
    r=r+round(1000/fs);
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function tQuantUnif

FPDF = PDFFn('Laplace');

Sym = true;
%Sym = false;
NlevA = [4 8 16 32 64 128 256];
for (Nlev = NlevA)
  [Yq, Xq, MSE, Entropy, SNRdB, sdV] = QuantUnif(Nlev, FPDF, Sym);
%  Xq
%  Yq
  fprintf('Nlev:%4d, SNR =%5.1f dB, s/V =%6.3f, Entropy =%5.2f bits\n', ...
      Nlev, SNRdB, sdV, Entropy);
end

return


text:
% test reading mouse data

%% data locations

% path to root directory containing the input data
inmoviedir = 'Y:\popcage_enriched\mousetrack_18';
intrxdir = 'Y:\adam\mousetrack_18\Results\Tracks';
% base name of the input experiment
inmovieexpname = 'b6_popcage_18_09.15.11_10.56.24.135';
intrxexpname = 'b6_popcage_18_2011.09.15_10.56.24.135';

% path to root directory to hold the output data
outdatadir = 'C:\Workspace\data';

%% parameters

% scale
mperpx = 0.00098387;
pxpermm = 1/(1000*mperpx);

% i assigned this arbitrarily for now
sex = {'F','M','M','F'};

% whether to make links or to copy
makelinks = true;

% frame interval length
frameintervallength = 108000;

% read in length of movie
inseqfile = fullfile(inmoviedir,[inmovieexpname,'.seq']);
headerinfo = r_readseqinfo(inseqfile);
nframes = headerinfo.m_iNumFrames;

%% main function call

for intervalstart = 1:frameintervallength:nframes,
  intervalend = intervalstart + frameintervallength - 1;
  if intervalend > nframes,
    break;
  end
  frameinterval = [intervalstart,intervalend];
  fprintf('Creating experiment for frames %d to %d...\n',intervalstart,intervalend);
  ConvertMouseTrx2Trx(inmoviedir,intrxdir,outdatadir,inmovieexpname,intrxexpname,...
    pxpermm,sex,'frameinterval',frameinterval,'makelinks',makelinks);
end

author: kristinbranson

repo: JAABA


text:
function done = save_model_expansion(method,goal_vexity,F_graph,properties)
global ALREADY_MODELLED REMOVE_THESE_IN_THE_END ALREADY_MODELLED_INDEX

done = 0;

for variables = properties.models
    save_place = find(variables == ALREADY_MODELLED_INDEX);
    if ~isempty(save_place)%length(ALREADY_MODELLED)>=variables & ~isempty(ALREADY_MODELLED{variables})
        % Ouch, we have already modelled this-one
        if strcmpi(ALREADY_MODELLED{save_place}.method,method)
            % Ok, already modelled using same approach
            done = 1;
            return
        elseif strcmpi(ALREADY_MODELLED{save_place}.method,'graph') & (strcmpi(method,'callback') | strcmpi(method,'integer') | strcmpi(method,'exact'))
            % Replace old graph model with exact model
            REMOVE_THESE_IN_THE_END = [REMOVE_THESE_IN_THE_END  ALREADY_MODELLED{save_place}.index];
            ALREADY_MODELLED{save_place}.goal_vexity = goal_vexity;
            ALREADY_MODELLED{save_place}.method = method;
            ALREADY_MODELLED{save_place}.index  = getlmiid(F_graph);
            ALREADY_MODELLED{save_place}.properties  = properties;
        elseif ((strcmpi(ALREADY_MODELLED{save_place}.method,'integer') | strcmpi(ALREADY_MODELLED{save_place}.method,'exact') | strcmpi(ALREADY_MODELLED{save_place}.method,'callback')) & strcmpi(method,'graph'))
            % Keep old stuff, we are done
            done = 1;
            return
        end
    else
        ALREADY_MODELLED{end+1}.goal_vexity = goal_vexity;
        ALREADY_MODELLED{end}.method = method;
        ALREADY_MODELLED{end}.index  = getlmiid(F_graph);
        ALREADY_MODELLED{end}.properties = properties;
        ALREADY_MODELLED_INDEX = [ALREADY_MODELLED_INDEX;variables];
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%  Figure 10.51      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
%  fig10_51.m is a script to generate Fig. 10.51, the rootlocus
%  for the fuel/air ratio with a linear sensor.  A (3,3)
%  Pade approximant is used  to approximate the delay.
clf;
% construct the F/A dynamics with the sensor time-constant too
f =[-50 ,    0,     0;
     0 ,   -1 ,    0;
    10 ,   10 ,  -10];

g =[25.0000;
    0.5000;
         0];
h =[0,     0,     1];
j=0;
[n3,d3]=pade(0.2,3); %  the delay Pade model

[np,dp]=ss2tf(f,g,h,j); % convert to polynomial form

np=[np(3:4)];  % remove the extraneous leading zeros

np=conv(np,n3);  % the plant numerator
dp =conv(dp,d3);  % the plant denominator
nc=[1, .3];
dc=[1, 0]; % the PI controller in polynomial form
% form the open-loop system
nol=conv(np,nc);
dol=conv(dp,dc);

rlocus(nol,dol);
v=[-30 2 -12 12];  % set the axes
axis(v);
grid;
title('Fig. 10.51  Rootlocus for the F/A control')


text:
function p = predictOneVsAll(all_theta, X)
%PREDICT Predict the label for a trained one-vs-all classifier. The labels
%are in the range 1..K, where K = size(all_theta, 1).
%  p = PREDICTONEVSALL(all_theta, X) will return a vector of predictions
%  for each example in the matrix X. Note that X contains the examples in
%  rows. all_theta is a matrix where the i-th row is a trained logistic
%  regression theta vector for the i-th class. You should set p to a vector
%  of values from 1..K (e.g., p = [1; 3; 1; 2] predicts classes 1, 3, 1, 2
%  for 4 examples)

m = size(X, 1);
num_labels = size(all_theta, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters (one-vs-all).
%               You should set p to a vector of predictions (from 1 to
%               num_labels).
%
% Hint: This code can be done all vectorized using the max function.
%       In particular, the max function can also return the index of the
%       max element, for more information see 'help max'. If your examples
%       are in rows, then, you can use max(A, [], 2) to obtain the max
%       for each row.
%

% =========================================================================

[probability indices] = max(sigmoid(all_theta * X'));
p = indices';

end

author: AvaisP


text:
%% UNO Rosenbrock
clc
%Objective
obj = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
%Setup Options
opts = optiset('solver','ipopt','display','iter');
%Build & Solve
Opt = opti('obj',obj,'ndec',2,'options',opts)
x0 = [0 0]';
[x,fval,exitflag,info]= solve(Opt,x0)
%Plot
% plot(Opt,[],1)

%% NLP1 Hock & Schittkowski #71
clc
%Objective & Gradient
obj = @(x) x(1)*x(4)*sum(x(1:3)) + x(3);
grad = @(x) [ x(1)*x(4) + x(4)*sum(x(1:3));
              x(1)*x(4);
              x(1)*x(4) + 1;
              x(1)*sum(x(1:3)) ];
%Linear Constraints
lb = ones(4,1);
ub = 5*ones(4,1);
%Nonlinear Constraints
nlcon = @(x) [ prod(x);
               sum(x.^2)];
nljac = @(x) [ prod(x)./x';
                2*x' ];
nlrhs = [25 40]';
nle = [1 0]'; % (>=, ==)
%Setup Options
opts = optiset('solver','ipopt','warnings','on','display','iter','solverOpts',ipoptset('linear_solver','pardiso'));
%Build & Solve
Opt = opti('obj',obj,'grad',grad,'nlmix',nlcon,nlrhs,nle,'nljac',nljac,'bounds',lb,ub,'options',opts)
x0 = [1 5 5 1]';
[x,fval,exitflag,info]= solve(Opt,x0)

info.Lambda

%%
clc
prob = amplRead('opf_662bus');
opts = optiset('solver','ipopt','solverOpts',ipoptset('linear_solver','pardiso'),'display','off');
opts2 = optiset('solver','ipopt','solverOpts',ipoptset('linear_solver','ma57'),'display','off');
opts3 = optiset('solver','ipopt','solverOpts',ipoptset('linear_solver','mumps'),'display','off');
O = opti(prob,opts)
O2 = opti(prob,opts2);
O3 = opti(prob,opts3);

[x,f,e,i] = solve(O);
[x2,f,e,i2] = solve(O2);
[x3,f,e,i3] = solve(O3);

norm(x-x2)

i
i2
i3

%%
clc
clear funcs opts

x0 = [0; 0];

funcs.objective = @(x) (1-x(1))^2 + 100 *(x(2)-x(1)^2)^2;
funcs.gradient = @(x)[2*x(1)-400*x(1)*(x(2)-x(1)^2)-2,200*x(2)-200*x(1)^2];

opts.ipopt.print_level = 5;
opts.ipopt.hessian_approximation = 'limited-memory';
opts.ipopt.ma57_pivot_order = 2;
opts.ipopt.linear_solver = 'pardiso';

[x,f] = ipopt(x0,funcs,opts)

f.eval

%%

prob = amplRead('ch3.nl')
opts = optiset('solver','ipopt','display','iter','solverOpts',ipoptset('linear_solver','pardiso'));

Opt = opti(prob,opts)

x = solve(Opt)

asl('close')

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function compute_mpii_crop_size(expidx)

% mpii 654
p = exp_params(expidx);

pidxs = p.pidxs;
num_joints = length(pidxs);
parts = get_parts();

% load annolist
load(p.trainGT);

num_images = length(annolist);

data = zeros(num_images, 6);

for i = 1:num_images
    if mod(i, 100) == 0
        fprintf('processing image %d/%d \n', i, num_images);
    end
    %filename = annolist(i).image.name;

    joint_list = get_anno_joints(annolist(i).annorect(1), p.pidxs, parts);
    joints = zeros(num_joints, 2);
    n = 0;
    for j = 1:num_joints
        jnt = joint_list(j, :);
        if ~isnan(jnt(1))
            n = n + 1;
            joints(n, :) = [jnt];
        end
    end
    joints = joints(1:n, :);
    data(i, :) = [mean(joints, 1) min(joints, [], 1) max(joints, [], 1)];
end

% [mean_x mean_y left_x top_y right_x bottom_y]

out_dir = fullfile(p.expDir, p.shortName, 'data', 'crop_sizes.mat');
save(out_dir, 'data');

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
%%
% Interpolation of 1D Gaussians

rep = 'results/gaussian-1d/';
[~,~] = mkdir(rep);

N = 1024;
x = linspace(0,1,N)';
normalize = @(x)x/sum(x(:));
gauss = @(m,s)normalize( exp( -(x-m).^2/(2*s^2) ) );

m0 = .1;  s0 = .025;
m1 = .8;  s1 = .08;

q = 16;
clf; hold on;
for i=1:q
    t = (i-1)/(q-1);
    m = (1-t)*m0+t*m1;
    s = (1-t)*s0+t*s1;
    lw = 2;
    plot(x, gauss(m,s), 'color', [1-t 0 t], 'LineWidth', lw);
end
set(gca, 'XTick', [], 'YTick', []);
box on; axis tight;
saveas(gcf, [rep 'interp-density.eps']);

clf; hold on;
for i=1:q
    t = (i-1)/(q-1);
    m = (1-t)*m0+t*m1;
    s = (1-t)*s0+t*s1;
    lw = 2;
    plot(x, gauss(m,s), 'color', [1-t 0 t], 'LineWidth', lw);
end
set(gca, 'XTick', [], 'YTick', []);
box on; axis tight;
saveas(gcf, [rep 'interp-density.eps'], 'epsc');

author: optimaltransport


text:
function plot_spreads(treedata,fig,lw,rel)

% Plots the spreads as a polar plot with different height layers presented
% with different colors. Inputs "fig" and "lw" define the figure number and
% the line width. Input Rel = 1 specifies relative spreads, i.e. the
% maximum spread is one, otherwise use the actual values.

if nargin == 2
    lw = 1;
    rel = 1;
elseif nargin == 3
    rel = 1;
end

spreads = treedata.spreads;
figure(fig)
n = size(spreads,1);
col = zeros(n,3);
col(:,1) = (0:1/n:(n-1)/n)';
col(:,3) = (1:-1/n:1/n)';
d = max(max(spreads));
D = [spreads(1,end) spreads(1,:)];
if rel
    polarplot(D/d,'-','Color',col(1,:),'Linewidth',lw)
else
    polarplot(D,'-','Color',col(1,:),'Linewidth',lw)
end
hold on
for i = 1:n
    D = [spreads(i,end) spreads(i,:)];
    if rel
        polarplot(D/d,'-','Color',col(i,:),'Linewidth',lw)
    else
        polarplot(D,'-','Color',col(i,:),'Linewidth',lw)
    end
end
hold off
if rel
    rlim([0 1])
else
    rlim([0 d])
end

author: InverseTampere

repo: TreeQSM


text:
function nii=nii2jnii(filename, format, varargin)
%
%    nii=nii2jnii(niifile, format, options)
%       or
%    nii2jnii(niifile, jniifile, options)
%    nii=nii2jnii(niifile)
%
%    A fast and portable NIFTI-1/2 and Analyze7.5 file parser and converter
%    to the text and binary JNIfTI formats defined in JNIfTI specification:
%    https://github.com/fangq/jnifti
%
%    This function is compatible with both MATLAB and GNU Octave.
%    It accepts .nii, .nii.gz, .hdr/.img and .hdr.gz/.img.gz input files
%
%    author: Qianqian Fang (q.fang <at> neu.edu)
%
%    input:
%        fname: the file name to the .nii, .nii.gz, .hdr/.img or .hdr.gz/.img.gz file
%        format:'nii' for reading the NIfTI-1/2/Analyze files;
%               'jnii' to convert the nii data into an in-memory JNIfTI structure.
%               'niiheader' return only the nii header without the image data
%
%               if format is not listed above and nii2jnii is called without
%               an output, format must be a string specifying the output JNIfTI
%               file name - *.jnii for text-based JNIfTI, or *.bnii for the
%               binary version
%        options: (optional) if saving to a .bnii file, please see the options for
%               savebj.m (part of JSONLab); if saving to .jnii, please see the
%               supported options for savejson.m (part of JSONLab).
%
%    output:
%        if the output is a JNIfTI data structure, it has the following subfield:
%          nii.NIFTIHeader -  a structure containing the 1-to-1 mapped NIFTI-1/2 header
%          nii.NIFTIData - the main image data array
%          nii.NIFTIExtension - a cell array contaiing the extension data buffers
%
%        when calling as nii=nii2jnii(file,'nii'), the output is a NIFTI object containing
%          nii.img: the data volume read from the nii file
%          nii.datatype: the data type of the voxel, in matlab data type string
%          nii.datalen: data count per voxel - for example RGB data has 3x
%                    uint8 per voxel, so datatype='uint8', datalen=3
%          nii.voxelbyte: total number of bytes per voxel: for RGB data,
%                    voxelbyte=3; also voxelbyte=header.bitpix/8
%          nii.hdr: file header info, a structure has the full nii header
%                    key subfileds include
%
%              sizeof_hdr: must be 348 (for NIFTI-1) or 540 (for NIFTI-2)
%              dim: short array, dim(2: dim(1)+1) defines the array size
%              datatype: the type of data stored in each voxel
%              bitpix: total bits per voxel
%              magic: must be 'ni1\0' or 'n+1\0'
%
%              For the detailed nii header, please see
%              https://nifti.nimh.nih.gov/pub/dist/src/niftilib/nifti1.h
%
%    dependency:
%          No external dependency if reading .nii/.hdr/.img files;
%
%          To load gzipped input files (.nii.gz/.hdr.gz/.img.gz), one must
%          install the ZMat Toolbox (http://github.com/fangq/zmat) and
%          JSONLab Toolbox (http://github.com/fangq/jsonlab);
%
%          To save files into the text/binary JNIfTI formatted files, one
%          need to install JSONLab (http://github.com/fangq/jsonlab).
%
%    this file was initially developed for the MCX project: https://github.com/fangq/mcx/blob/master/utils/mcxloadnii.m
%
%    this file is part of JNIfTI specification: https://github.com/fangq/jnifti
%
%    License: Apache 2.0, see https://github.com/fangq/jnifti for details
%

hdrfile=filename;
isnii=-1;
if(regexp(filename,'(\.[Hh][Dd][Rr](\.[Gg][Zz])*$|\.[Ii][Mm][Gg](\.[Gg][Zz])*$)'))
    isnii=0;
elseif(regexp(filename,'\.[Nn][Ii][Ii](\.[Gg][Zz])*$'))
    isnii=1;
end

if(isnii<0)
    error('file must be a NIfTI (.nii/.nii.gz) or Analyze 7.5 (.hdr/.img,.hdr.gz/.img.gz) data file');
end

if(regexp(filename,'\.[Ii][Mm][Gg](\.[Gg][Zz])*$'))
    hdrfile=regexprep(filename,'\.[Ii][Mm][Gg](\.[Gg][Zz])*$','.hdr$1');
end

niftiheader=niiformat('nifti1');

if(~isempty(regexp(hdrfile,'\.[Gg][Zz]$', 'once')) || (exist('OCTAVE_VERSION','builtin')~=0))
    finput=fopen(hdrfile,'rb');
    input=fread(finput,inf,'uint8=>uint8');
    fclose(finput);

    if(regexp(hdrfile,'\.[Gg][Zz]$'))
        if(~exist('gzipdecode','file'))
            error('To process zipped files, you must install gzipdecode.m from the JSONLab toolbox: http://github.com/fangq/jsonlab');
        end
        gzdata=gzipdecode(input);
    else
        gzdata=input;
    end
    clear input;
    nii.hdr=memmapstream(gzdata,niftiheader);
else
    fileinfo=dir(hdrfile);
    if(isempty(fileinfo))
        error('specified file does not exist');
    end
    header = memmapfile(hdrfile,             ...
       'Offset', 0,                          ...
       'Writable', false,                    ...
       'Format', niftiheader(1:end-(fileinfo.bytes<352),:));

    nii.hdr=header.Data(1);
end

[os,maxelem,dataendian]=computer;

if(nii.hdr.sizeof_hdr~=348 && nii.hdr.sizeof_hdr~=540)
    nii.hdr.sizeof_hdr=swapbytes(nii.hdr.sizeof_hdr);
end

if(nii.hdr.sizeof_hdr==540) % NIFTI-2 format
    niftiheader=niiformat('nifti2');
    if(exist('gzdata','var'))
        nii.hdr=memmapstream(gzdata,niftiheader);
    else
        header = memmapfile(hdrfile,                ...
         'Offset', 0,                           ...
         'Writable', false,                     ...
         'Format', niftiheader(1:end-(fileinfo.bytes<352),:));

         nii.hdr=header.Data(1);
    end
end

if(nii.hdr.dim(1)>7)
    names=fieldnames(nii.hdr);
    for i=1:length(names)
        nii.hdr.(names{i})=swapbytes(nii.hdr.(names{i}));
    end
    if(nii.hdr.sizeof_hdr>540)
        nii.hdr.sizeof_hdr=swapbytes(nii.hdr.sizeof_hdr);
    end
    if(dataendian=='B')
        dataendian='little';
    else
        dataendian='big';
    end
end

type2byte=[
        0  0  % unknown                      %
        1  0  % binary (1 bit/voxel)         %
        2  1  % unsigned char (8 bits/voxel) %
        4  2  % signed short (16 bits/voxel) %
        8  4  % signed int (32 bits/voxel)   %
       16  4  % float (32 bits/voxel)        %
       32  8  % complex (64 bits/voxel)      %
       64  8  % double (64 bits/voxel)       %
      128  3  % RGB triple (24 bits/voxel)   %
      255  0  % not very useful (?)          %
      256  1  % signed char (8 bits)         %
      512  2  % unsigned short (16 bits)     %
      768  4  % unsigned int (32 bits)       %
     1024  8  % long long (64 bits)          %
     1280  8  % unsigned long long (64 bits) %
     1536 16  % long double (128 bits)       %
     1792 16  % double pair (128 bits)       %
     2048 32  % long double pair (256 bits)  %
     2304  4  % 4 byte RGBA (32 bits/voxel)  %
];

type2str={
    'uint8'    0   % unknown                       %
    'uint8'    0   % binary (1 bit/voxel)          %
    'uint8'    1   % unsigned char (8 bits/voxel)  %
    'uint16'   1   % signed short (16 bits/voxel)  %
    'int32'    1   % signed int (32 bits/voxel)    %
    'single'   1   % float (32 bits/voxel)         %
    'single'   2   % complex (64 bits/voxel)       %
    'double'   1   % double (64 bits/voxel)        %
    'uint8'    3   % RGB triple (24 bits/voxel)    %
    'uint8'    0   % not very useful (?)           %
    'int8'     1   % signed char (8 bits)          %
    'uint16'   1   % unsigned short (16 bits)      %
    'uint32'   1   % unsigned int (32 bits)        %
    'int64'    1   % long long (64 bits)           %
    'uint64'   1   % unsigned long long (64 bits)  %
    'uint8'    16  % long double (128 bits)        %
    'uint8'    16  % double pair (128 bits)        %
    'uint8'    32  % long double pair (256 bits)   %
    'uint8'    4   % 4 byte RGBA (32 bits/voxel)   %
};

typeidx=find(type2byte(:,1)==nii.hdr.datatype);

nii.datatype=type2str{typeidx,1};
nii.datalen=type2str{typeidx,2};
nii.voxelbyte=type2byte(typeidx,2);

if(type2byte(typeidx,2)==0)
    nii.img=[];
    return;
end

if(type2str{typeidx,2}>1)
    nii.hdr.dim=[nii.hdr.dim(1)+1 uint16(nii.datalen) nii.hdr.dim(2:end)];
end

if(nargin>1 && strcmp(format,'niiheader'))
    return;
end

if(regexp(filename,'\.[Hh][Dd][Rr](\.[Gg][Zz])*$'))
    filename=regexprep(filename,'\.[Hh][Dd][Rr](\.[Gg][Zz])*$','.img$1');
end

imgbytenum=prod(nii.hdr.dim(2:nii.hdr.dim(1)+1))*nii.voxelbyte;

if(isnii==0 && ~isempty(regexp(filename,'\.[Gg][Zz]$', 'once')))
    finput=fopen(filename,'rb');
    input=fread(finput,inf,'uint8=>uint8');
    fclose(finput);
    gzdata=gzipdecode(input);
    nii.img=typecast(gzdata(1:imgbytenum),nii.datatype);
else
    if(~exist('gzdata','var'))
        fid=fopen(filename,'rb');
        if(isnii)
            fseek(fid,nii.hdr.vox_offset,'bof');
        end
        nii.img=fread(fid,imgbytenum,[nii.datatype '=>' nii.datatype]);
        fclose(fid);
    else
        nii.img=typecast(gzdata(nii.hdr.vox_offset+1:nii.hdr.vox_offset+imgbytenum),nii.datatype);
    end
end

nii.img=reshape(nii.img,nii.hdr.dim(2:nii.hdr.dim(1)+1));

if(nargin>1 && strcmp(format,'nii'))
    return;
end

nii0=nii;
nii=struct();
nii.NIFTIHeader.NIIHeaderSize=  nii0.hdr.sizeof_hdr;
if(isfield(nii0.hdr,'data_type'))
    nii.NIFTIHeader.A75DataTypeName=deblank(char(nii0.hdr.data_type));
    nii.NIFTIHeader.A75DBName=      deblank(char(nii0.hdr.db_name));
    nii.NIFTIHeader.A75Extends=     nii0.hdr.extents;
    nii.NIFTIHeader.A75SessionError=nii0.hdr.session_error;
    nii.NIFTIHeader.A75Regular=     nii0.hdr.regular;
end
nii.NIFTIHeader.DimInfo.Freq=   bitand(nii0.hdr.dim_info,7);
nii.NIFTIHeader.DimInfo.Phase=  bitand(bitshift(nii0.hdr.dim_info,-3),7);
nii.NIFTIHeader.DimInfo.Slice=  bitand(bitshift(nii0.hdr.dim_info,-6),7);
nii.NIFTIHeader.Dim=            nii0.hdr.dim(2:2+nii0.hdr.dim(1)-1);
nii.NIFTIHeader.Param1=         nii0.hdr.intent_p1;
nii.NIFTIHeader.Param2=         nii0.hdr.intent_p2;
nii.NIFTIHeader.Param3=         nii0.hdr.intent_p3;
nii.NIFTIHeader.Intent=         niicodemap('intent',nii0.hdr.intent_code);
nii.NIFTIHeader.DataType=       niicodemap('datatype',nii0.hdr.datatype);
nii.NIFTIHeader.BitDepth=       nii0.hdr.bitpix;
nii.NIFTIHeader.FirstSliceID=   nii0.hdr.slice_start;
nii.NIFTIHeader.VoxelSize=      nii0.hdr.pixdim(2:2+nii0.hdr.dim(1)-1);
nii.NIFTIHeader.Orientation=    struct('x','r','y','a','z','s');
if(nii0.hdr.pixdim(1)<0)
    nii.NIFTIHeader.Orientation=    struct('x','l','y','a','z','s');
end
nii.NIFTIHeader.NIIByteOffset=  nii0.hdr.vox_offset;
nii.NIFTIHeader.ScaleSlope=     nii0.hdr.scl_slope;
nii.NIFTIHeader.ScaleOffset=    nii0.hdr.scl_inter;
nii.NIFTIHeader.LastSliceID=    nii0.hdr.slice_end;
nii.NIFTIHeader.SliceType=      niicodemap('slicetype',nii0.hdr.slice_code);
nii.NIFTIHeader.Unit.L=         niicodemap('unit',bitand(nii0.hdr.xyzt_units, 7));
nii.NIFTIHeader.Unit.T=         niicodemap('unit',bitand(nii0.hdr.xyzt_units, 56));
nii.NIFTIHeader.MaxIntensity=   nii0.hdr.cal_max;
nii.NIFTIHeader.MinIntensity=   nii0.hdr.cal_min;
nii.NIFTIHeader.SliceTime=      nii0.hdr.slice_duration;
nii.NIFTIHeader.TimeOffset=     nii0.hdr.toffset;
if(isfield(nii0.hdr,'glmax'))
    nii.NIFTIHeader.A75GlobalMax=   nii0.hdr.glmax;
    nii.NIFTIHeader.A75GlobalMin=   nii0.hdr.glmin;
end
nii.NIFTIHeader.Description=    deblank(char(nii0.hdr.descrip));
nii.NIFTIHeader.AuxFile=        deblank(char(nii0.hdr.aux_file));
nii.NIFTIHeader.QForm=          nii0.hdr.qform_code;
nii.NIFTIHeader.SForm=          nii0.hdr.sform_code;
nii.NIFTIHeader.Quatern.b=      nii0.hdr.quatern_b;
nii.NIFTIHeader.Quatern.c=      nii0.hdr.quatern_c;
nii.NIFTIHeader.Quatern.d=      nii0.hdr.quatern_d;
nii.NIFTIHeader.QuaternOffset.x=nii0.hdr.qoffset_x;
nii.NIFTIHeader.QuaternOffset.y=nii0.hdr.qoffset_y;
nii.NIFTIHeader.QuaternOffset.z=nii0.hdr.qoffset_z;
nii.NIFTIHeader.Affine(1,:)=    nii0.hdr.srow_x;
nii.NIFTIHeader.Affine(2,:)=    nii0.hdr.srow_y;
nii.NIFTIHeader.Affine(3,:)=    nii0.hdr.srow_z;
nii.NIFTIHeader.Name=           deblank(char(nii0.hdr.intent_name));
nii.NIFTIHeader.NIIFormat=      deblank(char(nii0.hdr.magic));
if(isfield(nii0.hdr,'extension'))
    nii.NIFTIHeader.NIIExtender=    nii0.hdr.extension;
end
nii.NIFTIHeader.NIIQfac_=       nii0.hdr.pixdim(1);
nii.NIFTIHeader.NIIEndian_=     dataendian;
if(isfield(nii0.hdr,'reserved'))
    nii.NIFTIHeader.NIIUnused_= nii0.hdr.reserved;
end

nii.NIFTIData=nii0.img;

if(isfield(nii0.hdr,'extension') && nii0.hdr.extension(1)>0)
    fid=fopen(filename,'rb');
    fseek(fid,nii0.hdr.sizeof_hdr+4,'bof');
    nii.NIFTIExtension=cell(1);
    count=1;
    while(ftell(fid)<nii0.hdr.vox_offset)
       nii.NIFTIExtension{count}.Size=fread(fid,1,'int32=>int32');
       nii.NIFTIExtension{count}.Type=fread(fid,1,'int32=>int32');
       if(strcmp(dataendian,'big'))
           nii.NIFTIExtension{count}.Size=swapbytes(nii.NIFTIExtension{count}.Size);
           nii.NIFTIExtension{count}.Type=swapbytes(nii.NIFTIExtension{count}.Type);
       end
       nii.NIFTIExtension{count}.x0x5F_ByteStream_=fread(fid,nii.NIFTIExtension{count}.Size-8,'uint8=>uint8');
       count=count+1;
    end
    fclose(fid);
end

if(nargout==0 && strcmp(format,'nii')==0 && strcmp(format,'jnii')==0)
    if(~exist('savejson','file'))
        error('you must first install JSONLab from http://github.com/fangq/jsonlab/');
    end
    if(regexp(format,'\.jnii$'))
        savejson('',nii,'FileName',format,varargin{:});
    elseif(regexp(format,'\.bnii$'))
        savebj('',nii,'FileName',format,varargin{:});
    else
        error('file suffix must be .jnii for text JNIfTI or .bnii for binary JNIfTI');
    end
end


text:
%
% step 3: evaluation of performance
%
% using our third set of prices, we estimate dp at each time interval,
% if dp > t and current position <= 0 , we buy
% if dp < -t and current position >= 0, we sell
% else, do:nothing
%

% trade using the above algorithm. returns expected profit
%
% given a list of prices
% assumes k-means clustered patterns have already been calculated
% assumes parameters w_i have already been calculated
% position is 0 or 1 (we have nothing, we have a bitcoin)
% bank is the amount of cash we have
% threshold is the threshold for buying/selling
% defined in the paper
function [error,jinzhi,bank,buy,sell,proba] = brtrade(prices, bidVolume,askVolume, fee)
    assert(exist('thetas.mat','file')==2)
    load('thetas.mat');
    assert(isequal(length(prices), length(bidVolume)));
    assert(isequal(length(prices), length(askVolume)));
    position = 0;
    bank = 0;
    jinzhi = zeros(length(prices)-750, 1);
    error = 0;
    %current error metric is sum(abs(error))/time interval = ~.9
    %current error = 0.06
    buy = [];
    sell = [];
    counttotal = 0;
    counts = 0;
    temp = 0;
    for t = 720:length(prices)-1
        price180 = zscore(prices(t-179:t));
        price360 = zscore(prices(t-359:t));
        price720 = zscore(prices(t-719:t));

        %average price change dp_j is given by bayesian regression
        dp1 = bayesian(price180,kmeans180s);
        dp2 = bayesian(price360,kmeans360s);
        dp3 = bayesian(price720,kmeans720s);

        r = (bidVolume(t)-askVolume(t))/(bidVolume(t)+askVolume(t));

        dp = theta0 +  theta(1)*dp1 + theta(2)*dp2 + theta(3)*dp3;% + theta(4)*r;

        % compare price at t+1 with predicted price jump
        error = error + abs(prices(t+1)-prices(t)-dp);

        % calculate transaction fee??
        % threshold 1 and 2 before...but
        % there is definitely not going to be a
        % 5-8$ price jump predicted in the next ten seconds
        % need to consider TODO
		fee = 0;
		if (fee == 0)
            tfee_buy = 0.001;
            tfee_sell = 0.003;
        else
            tfee_buy = fee*prices(2)/100;
            tfee_sell = tfee_buy;
        end
        %BUY
        if (dp > tfee_buy && position == 0)
            position = 1;
            temp = prices(t);
            fprintf('Buying at %d\n', temp);
            buy = [buy;t];
        end
        %SELL
        if (dp < -tfee_sell && position == 1)
            position = 0;
            bank = bank + prices(t)-temp;
            fprintf('Selling at %d\n', prices(t));
            sell = [sell;t];
            counttotal = counttotal+1;
            if prices(t)-temp>0
                counts = counts+1;
            end
        end

        jinzhi(t) = bank;
    end

    % forces us to close the position at the end
    % tradeoffs to this decision
    % on one side: more realistic
    % but the algorithm doesn't yet account for it
    if (position == 1)
        bank = bank + prices(t)-temp;
        fprintf('Final sale at %d\n', prices(t));
        sell = [sell;t];
        counttotal = counttotal+1;
        if prices(t)-temp>0
            counts = counts+1;
        end
    end
    proba = (counts./counttotal)*100;
    end

author: panditanvita

repo: BTCpredictor


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

%*****************************************************************************80
%
%% EMPIRICAL_DISCRETE_PDF evaluates the Empirical Discrete PDF.
%
%  Discussion:
%
%    A set of A values C(1:A) are assigned nonnegative weights B(1:A),
%    with at least one B nonzero.  The probability of C(I) is the
%    value of B(I) divided by the sum of the weights.
%
%    The C's must be distinct, and given in ascending order.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%
%    Input, integer A, the number of values.
%    0 < A.
%
%    Input, real B(A), the weights of each value.
%    0 <= B(1:A) and at least one value is nonzero.
%
%    Input, real C(A), the values.
%    The values must be distinct and in ascending order.
%
%    Output, real PDF, the value of the PDF.
%
  for i = 1 : a
    if ( x == c(i) )
      pdf = b(i) / sum ( b(1:a) );
      return
    end
  end

  pdf = 0.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Average European Diet defintion. For details, please see https://www.vmh.life/#nutrition
% Ines Thiele 2016-2019

Diet = {%'Exchange reactions'	'Total in mmol/(human*day)'
%'reaction'	'fluxValue'
'EX_10fthf(e)'	'0.000125752'
'EX_5mthf(e)'	'0.000129311'
'EX_adpcbl(e)'	'3.61E-06'
'EX_ala_L(e)'	'37.20856362'
'EX_arach(e)'	'1.033055397'
'EX_arachd(e)'	'0.536810751'
'EX_arg_L(e)'	'20.3518771'
'EX_ascb_L(e)'	'0.583441166'
'EX_asp_L(e)'	'47.08674913'
'EX_avite1(e)'	'0.03685808'
'EX_btn(e)'	'0.000178666'
'EX_but(e)'	'26.33425644'
'EX_ca2(e)'	'18.5490793'
'EX_caro(e)'	'0.003585581'
'EX_CE2510(e)'	'1.086665034'
'EX_CE4843(e)'	'0.114865263'
'EX_cellul(e)'	'0.052346876'
'EX_chsterol(e)'	'1.23218851'
'EX_cl(e)'	'187.7787493'
'EX_clpnd(e)'	'0.048376886'
'EX_crvnc(e)'	'0.16630008'
'EX_cu2(e)'	'0.035674944'
'EX_cys_L(e)'	'7.511831638'
'EX_dca(e)'	'10.74842619'
'EX_ddca(e)'	'13.45904717'
'EX_doco13ac(e)'	'0.560730443'
'EX_docosac(e)'	'0.454096231'
'EX_etoh(e)'	'234.4340155'
'EX_fe2(e)'	'0.109248814'
'EX_fe3(e)'	'0.109248814'
'EX_fol(e)'	'0.000134614'
'EX_fru(e)'	'106.5058652'
'EX_gal(e)'	'0.029974039'
'EX_glc_D(e)'	'135.5803501'
'EX_glu_L(e)'	'103.1222005'
'EX_gly(e)'	'36.76468629'

'EX_h2o(e)'	'165892.3425'
'EX_hdca(e)'	'92.78506719'
'EX_hdcea(e)'	'7.120473022'
'EX_his_L(e)'	'14.02407663'
'EX_hpdca(e)'	'1.779524923'
'EX_i(e)'	'0.00201238'
'EX_ile_L(e)'	'28.92167513'
'EX_k(e)'	'83.43380658'
'EX_lcts(e)'	'4.126714705'
'EX_leu_L(e)'	'46.88780991'
'EX_lgnc(e)'	'0.222615856'
'EX_lnlc(e)'	'54.98546187'
'EX_lnlnca(e)'	'4.163336898'
'EX_lys_L(e)'	'33.31154825'
'EX_malt(e)'	'8.413758248'
'EX_met_L(e)'	'10.72740469'
'EX_mg2(e)'	'14.52993211'
'EX_mn2(e)'	'0.074702331'
'EX_mnl(e)'	'0.170169038'
'EX_na1(e)'	'172.4462717'
'EX_nac(e)'	'0.308616192'
'EX_ncam(e)'	'0.1455276'
'EX_ocdca(e)'	'35.85086639'
'EX_ocdcea(e)'	'101.5927326'
'EX_octa(e)'	'6.845572908'
'EX_phe_L(e)'	'20.71498664'
'EX_phyQ(e)'	'0.000136877'
'EX_pi(e)'	'14.22660928'
'EX_pnto_R(e)'	'0.025569705'
'EX_pro_L(e)'	'53.81527918'
'EX_ptdca(e)'	'2.89407663'
'EX_pydam(e)'	'0.00317177'
'EX_pydx(e)'	'0.003210459'
'EX_pydxn(e)'	'0.003172205'
'EX_retinol(e)'	'3.002252387'
'EX_ribflv(e)'	'0.003693234'
'EX_sbt_D(e)'	'0.032935943'
'EX_ser_L(e)'	'37.37922556'
'EX_strch1(e)'	'33.03699973'
'EX_strdnc(e)'	'0'
'EX_sucr(e)'	'202.3830217'
'EX_thf(e)'	'0.000133696'
'EX_thm(e)'	'5.407858013'
'EX_thr_L(e)'	'25.88105024'
'EX_trp_L(e)'	'4.554824772'
'EX_ttdca(e)'	'30.21138004'
'EX_ttdcea(e)'	'3.783984795'
'EX_tyr_L(e)'	'15.64685397'
'EX_urate(e)'	'2.919392803'
'EX_val_L(e)'	'37.78437731'
'EX_vitd3(e)'	'8.76E-06'
'EX_xylt(e)'	'0.141969085'
'EX_zn2(e)'	'0.150963597'
'EX_adrn(e)'	'0'

% add bile acids
'EX_gchola(e)'	'10'
'EX_tdchola(e)'	'10'
'EX_tchola(e)'	'10'
'EX_dgchol(e)'	'10'
'EX_M01989(e)'	'10'
'EX_M01989(e)'	'10'

    };
Diet = regexprep(Diet,'EX_','Diet_EX_');
Diet = regexprep(Diet,'\(e\)','\[d\]');

author: opencobra

repo: cobratoolbox


text:
function c=ref_edgt(f,g,a,M)
%REF_EDGT   Reference Even Discrete Gabor transform
%   Usage  c=ref_edgt(f,g,a,M);
%
%   The input window must be odd-centered.

L=size(f,1);
W=size(f,2);

N=L/a;
M=L/b;

F=zeros(L,M*N);

l=(0:L-1)';
for n=0:N-1
  for m=0:M-1
    F(:,1+m+n*M)=exp(2*pi*i*m.*(l+.5)*b/L).*circshift(g,n*a);
  end;
end;

c=F'*f;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function showCollectionStatistics(collections, legendString, ...
    collectionColors, collectionMarkers)
% Displays the items in collectionStats uusing scatter plots

%% Plot median versus SDR overall channel deviation
baseTitle = 'Median versus SDR overall channel deviation';
figure ('Name', baseTitle, 'Color', [1, 1, 1]);
hold on
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    plot(stats(:, s.medDev), stats(:, s.rSDDev), ...
        'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'LineStyle', 'none', ...
        'MarkerSize', 12, 'LineWidth', 2);
end
xlabel('Median overall channel deviation');
ylabel('Robust SD overall channel deviation');
title(baseTitle, 'interpreter', 'none')
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    medDev = median(stats(:, s.medDev));
    rSDDev = median(stats(:, s.rSDDev));
    plot(medDev, rSDDev, 'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'MarkerSize', 14, 'LineWidth', 3);
end
legend(legendString, 'Location', 'NorthEast')
box on
hold off

%% Plot median versus SDR window channel deviation
baseTitle = 'Median versus SDR window channel deviation';
figure ('Name', baseTitle, 'Color', [1, 1, 1]);
hold on
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    plot(stats(:, s.medWinDev), stats(:, s.rSDWinDev), ...
        'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'LineStyle', 'none', ...
        'MarkerSize', 12, 'LineWidth', 2);
end
xlabel('Median window channel deviation');
ylabel('Robust SD window channel deviation');
title(baseTitle, 'interpreter', 'none')
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    medDev = median(stats(:, s.medWinDev));
    rSDDev = median(stats(:, s.rSDWinDev));
    plot(medDev, rSDDev, 'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'MarkerSize', 14, 'LineWidth', 3);
end
legend(legendString, 'Location', 'NorthEast')
box on
hold off

%% Plot median maximum correlation versus median window channel deviation
baseTitle = 'Median max correlation versus median window channel deviation';
figure ('Name', baseTitle, 'Color', [1, 1, 1]);
hold on
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    plot(stats(:, s.medCor), stats(:, s.medWinDev), ...
        'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'LineStyle', 'none', ...
        'MarkerSize', 12, 'LineWidth', 2);
end
xlabel('Median maximum correlation');
ylabel('Median window channel deviation');
title(baseTitle, 'interpreter', 'none')
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    medDev = median(stats(:, s.medWinDev));
    medCor = median(stats(:, s.medCor));
    plot(medCor, medDev, 'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'MarkerSize', 14, 'LineWidth', 3);
end
legend(legendString, 'Location', 'NorthEast')
box on
hold off

%% Plot median maximum correlation versus window channel deviation ratio
baseTitle = 'Median max correlation versus window channel deviation ratio';
figure ('Name', baseTitle, 'Color', [1, 1, 1]);
hold on
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    devRatio = stats(:, s.rSDWinDev)./stats(:, s.medWinDev);
    plot(stats(:, s.medCor), devRatio, ...
        'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'LineStyle', 'none', ...
        'MarkerSize', 12, 'LineWidth', 2);
end
xlabel('Median maximum correlation');
ylabel('Window channel deviation ratio');
title(baseTitle, 'interpreter', 'none')
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    medRatio = median(stats(:, s.rSDWinDev)./stats(:, s.medWinDev));
    medCor = median(stats(:, s.medCor));
    plot(medCor, medRatio, 'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'MarkerSize', 14, 'LineWidth', 3);
end
legend(legendString, 'Location', 'NorthEast')
box on
hold off

%% Plot average maximum correlation versus median window channel deviation
baseTitle = 'Average max correlation versus median window channel deviation';
figure ('Name', baseTitle, 'Color', [1, 1, 1]);
hold on
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    plot(stats(:, s.aveCor), stats(:, s.medWinDev), ...
        'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'LineStyle', 'none', ...
        'MarkerSize', 12, 'LineWidth', 2);
end
xlabel('Average maximum correlation');
ylabel('Median window channel deviation');
title(baseTitle, 'interpreter', 'none')
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    medDev = median(stats(:, s.medWinDev));
    aveCor = median(stats(:, s.aveCor));
    plot(aveCor, medDev, 'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'MarkerSize', 14, 'LineWidth', 3);
end
legend(legendString, 'Location', 'NorthEast')
box on
hold off

%% Plot average maximum correlation versus window channel deviation ratio
baseTitle = 'Average max correlation versus window channel deviation ratio';
figure ('Name', baseTitle, 'Color', [1, 1, 1]);
hold on
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    devRatio = stats(:, s.rSDWinDev)./stats(:, s.medWinDev);
    plot(stats(:, s.aveCor), devRatio, ...
        'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'LineStyle', 'none', ...
        'MarkerSize', 12, 'LineWidth', 2);
end
xlabel('Average maximum correlation');
ylabel('Window channel deviation ratio');
title(baseTitle, 'interpreter', 'none')
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    medRatio = median(stats(:, s.rSDWinDev)./stats(:, s.medWinDev));
    aveCor = median(stats(:, s.aveCor));
    plot(aveCor, medRatio, 'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'MarkerSize', 14, 'LineWidth', 3);
end
legend(legendString, 'Location', 'NorthEast')
box on
hold off

%% Plot median versus SDR overall HF noise
baseTitle = 'Median versus SDR overall HF noise';
figure ('Name', baseTitle, 'Color', [1, 1, 1]);
hold on
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    plot(stats(:, s.medHF), stats(:, s.rSDHF), ...
        'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'LineStyle', 'none', ...
        'MarkerSize', 12, 'LineWidth', 2);
end
xlabel('Median overall HF noise');
ylabel('Robust SD overall HF noise');
title(baseTitle, 'interpreter', 'none')
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    medHF = median(stats(:, s.medHF));
    rSDHF = median(stats(:, s.rSDHF));
    plot(medHF, rSDHF, 'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'MarkerSize', 14, 'LineWidth', 3);
end
legend(legendString, 'Location', 'NorthEast')
box on
hold off

%% Plot median versus SDR window HF noise
baseTitle = 'Median versus SDR window HF noise';
figure ('Name', baseTitle, 'Color', [1, 1, 1]);
hold on
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    plot(stats(:, s.medWinHF), stats(:, s.rSDWinHF), ...
        'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'LineStyle', 'none', ...
        'MarkerSize', 12, 'LineWidth', 2);
end
xlabel('Median window HF noise');
ylabel('Robust SD window HF noise');
title(baseTitle, 'interpreter', 'none')
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    medHF = median(stats(:, s.medWinHF));
    rSDHF = median(stats(:, s.rSDWinHF));
    plot(medHF, rSDHF, 'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'MarkerSize', 14, 'LineWidth', 3);
end
legend(legendString, 'Location', 'NorthEast')
box on
hold off

%% Plot median maximum correlation versus median window HF noise
baseTitle = 'Median max correlation versus median window HF noise';
figure ('Name', baseTitle, 'Color', [1, 1, 1]);
hold on
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    plot(stats(:, s.medCor), stats(:, s.medWinHF), ...
        'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'LineStyle', 'none', ...
        'MarkerSize', 12, 'LineWidth', 2);
end
xlabel('Median maximum correlation');
ylabel('Median window channel HF noise');
title(baseTitle, 'interpreter', 'none')
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    medHF = median(stats(:, s.medWinHF));
    medCor = median(stats(:, s.medCor));
    plot(medCor, medHF, 'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'MarkerSize', 14, 'LineWidth', 3);
end
legend(legendString, 'Location', 'NorthEast')
box on
hold off

%% Plot median maximum correlation versus window channel HF noise ratio
baseTitle = 'Median max correlation versus window channel HF noise ratio';
figure ('Name', baseTitle, 'Color', [1, 1, 1]);
hold on
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    HFRatio = stats(:, s.rSDWinHF)./stats(:, s.medWinHF);
    plot(stats(:, s.medCor), HFRatio, ...
        'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'LineStyle', 'none', ...
        'MarkerSize', 12, 'LineWidth', 2);
end
xlabel('Median maximum correlation');
ylabel('Window channel HF noise ratio');
title(baseTitle, 'interpreter', 'none')
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    medRatio = median(stats(:, s.rSDWinHF)./stats(:, s.medWinHF));
    medCor = median(stats(:, s.medCor));
    plot(medCor, medRatio, 'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'MarkerSize', 14, 'LineWidth', 3);
end
legend(legendString, 'Location', 'NorthEast')
box on
hold off

%% Plot average maximum correlation versus median window channel HF noise
baseTitle = 'Average max correlation versus median window channel HF noise';
figure ('Name', baseTitle, 'Color', [1, 1, 1]);
hold on
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    plot(stats(:, s.aveCor), stats(:, s.medWinHF), ...
        'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'LineStyle', 'none', ...
        'MarkerSize', 12, 'LineWidth', 2);
end
xlabel('Average maximum correlation');
ylabel('Median window channel HF noise');
title(baseTitle, 'interpreter', 'none')
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    medHF = median(stats(:, s.medWinHF));
    aveCor = median(stats(:, s.aveCor));
    plot(aveCor, medHF, 'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'MarkerSize', 14, 'LineWidth', 3);
end
legend(legendString, 'Location', 'NorthEast')
box on
hold off

%% Plot average maximum correlation versus window channel HF noise ratio
baseTitle = 'Average max correlation versus window channel HF noise ratio';
figure ('Name', baseTitle, 'Color', [1, 1, 1]);
hold on
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    HFRatio = stats(:, s.rSDWinHF)./stats(:, s.medWinHF);
    plot(stats(:, s.aveCor), HFRatio, ...
        'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'LineStyle', 'none', ...
        'MarkerSize', 12, 'LineWidth', 2);
end
xlabel('Average maximum correlation');
ylabel('Window channel HF noise ratio');
title(baseTitle, 'interpreter', 'none')
for k = 1:length(collections)
    stats = collections{k}.statistics;
    s = collections{k}.statisticsIndex;
    medRatio = median(stats(:, s.rSDWinHF)./stats(:, s.medWinHF));
    aveCor = median(stats(:, s.aveCor));
    plot(aveCor, medRatio, 'Color', collectionColors(k, :), ...
        'Marker', collectionMarkers{k}, 'MarkerSize', 14, 'LineWidth', 3);
end
legend(legendString, 'Location', 'NorthEast')
box on
hold off

author: VisLab

repo: EEG-Clean-Tools


text:
 % Copyright 2001, Brown University, Providence, Rhode Island.
 %
 % All Rights Reserved
 %
 % Permission to use this software for noncommercial research and
 % educational purposes is hereby granted without fee.
 % Redistribution, sale, or incorporation of this software into a
 % commercial product is prohibited.
 %
 % BROWN UNIVERSITY DISCLAIMS ANY AND ALL WARRANTIES WITH REGARD TO
 % THIS SOFTWARE,INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY
 % AND FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL BROWN
 % UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
 % DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE,
 % DATA OR PROFITS.

  function [a,b] = am282ab(x,y)

   dims = size(x);
   Np = dims(1)*dims(2);

   a = zeros(dims);
   b = zeros(dims);
   for n=1:Np  % hack
    if(y(n) ~= 1)
     a(n) = 2*(1+x(n))/(1-y(n))-1;
    else
     a(n) = -1;
    end
   end

   b = y;

repo: Aero-matlab


text:
function ns = count_squares(I,x1,y1,x2,y2,win);

[ny,nx] = size(I);

if ((x1-win <= 0) || (x1+win >= nx) || (y1-win <= 0) || (y1+win >= ny) || ...
        (x2-win <= 0) || (x2+win >= nx) || (y2-win <= 0) || (y2+win >= ny))
    ns = -1;
    return;
end;

if ((x1 - x2)^2+(y1-y2)^2) <  win,
    ns = -1;
    return;
end;

lambda = [y1 - y2;x2 - x1;x1*y2 - x2*y1];
lambda = 1/sqrt(lambda(1)^2 + lambda(2)^2) * lambda;
l1 = lambda + [0;0;win];
l2 = lambda - [0;0;win];
dx = x2-x1;
dy = y2 - y1;

if abs(dx) > abs(dy),
   if x2 > x1,
      xs = x1:x2;
   else
      xs = x1:-1:x2;
   end;
   ys = -(lambda(3) + lambda(1)*xs)/lambda(2);
else
   if y2 > y1,
       ys = y1:y2;
   else
       ys = y1:-1:y2;
   end;
   xs = -(lambda(3) + lambda(2)*ys)/lambda(1);
end;

Np = length(xs);
xs_mat = ones(2*win + 1,1)*xs;
ys_mat = ones(2*win + 1,1)*ys;
win_mat = (-win:win)'*ones(1,Np);
xs_mat2 = round(xs_mat - win_mat * lambda(1));
ys_mat2 = round(ys_mat - win_mat * lambda(2));
ind_mat = (xs_mat2 - 1) * ny + ys_mat2;
ima_patch = zeros(2*win + 1,Np);
ima_patch(:) = I(ind_mat(:));

%ima2 = ima_patch(:,win+1:end-win);

filtk = [ones(win,Np);zeros(1,Np);-ones(win,Np)];
out_f = sum(filtk.*ima_patch);
out_f_f = conv2(out_f,[1/4 1/2 1/4],'same');
out_f_f = out_f_f(win+1:end-win);
ns = length(find(((out_f_f(2:end)>=0)&(out_f_f(1:end-1)<0)) | ((out_f_f(2:end)<=0)&(out_f_f(1:end-1)>0))))+1;

return;

repo: ekfmonoslam


text:
function y = Insert(x,lines,where)

%Insert - Insert lines in a matrix.
%
%  USAGE
%
%    result = Insert(matrix,lines,indices)
%
%    array          matrix where the lines should be inserted
%    lines          list of values to insert
%    indices        list of (possibly repeated) matrix line numbers
%                   after which new lines should be inserted (use 0
%                   to insert before first line)
%
%  EXAMPLES
%
%    >> Insert([1;2;3;4;5],[-1;-2;-3],[0;3;5])
%
%    ans =
%
%      -1
%       1
%       2
%       3
%      -2
%       4
%       5
%      -3

% Copyright (C) 2004-2011 by Michaël Zugaro
%
% 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.

where = where(:);
mx = size(x,1);
nx = size(x,2);
ml = length(where);
if size(lines,1) == 1,
	lines = repmat(lines,ml,1);
end

% Lines in x must be moved to new indices in y
w = Accumulate(where+1,1,[mx+1,1]);
new = cumsum(w)+(1:mx+1)';
new(end) = [];

% Copy x into y
y = zeros(mx+ml,nx);
y(new,:) = x;

inserted = where+(1:ml)';
y(inserted,:) = lines;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
function showcam(c, T, P)

    if nargin < 2
        T = c.T;
    end

    T = SE3(T);

    clf
    axis([-3, 3, -3, 3, -3 3])
    daspect([1 1 1])
    hold on
    grid
    [x,y,z] = sphere(20);

    hg = hgtransform;
    surf(x,y,z, 'FaceColor', [0.8 0.8 1], 'EdgeColor', 0.5*[0.8 0.8 1], ...
        'EdgeLighting', 'gouraud', 'Parent', hg)
    light
    lighting gouraud
    set(hg, 'Matrix', T.T);

    trplot(T, 'length', 1.6, 'arrow')

    axis
    limits = reshape(axis, 2, []);
    maxdim = max(diff(limits));

    o = T.t;

    if nargin > 2
        for i=1:numcols(P)
            plot3([o(1) P(1,i)], [o(2) P(2,i)], [o(3) P(3,i)], 'r');
            plot_sphere(P(:,i), maxdim*0.02, 'r');
        end
    end
    hold off
    xlabel('x'); ylabel('y'); zlabel('z');

author: petercorke

repo: machinevision-toolbox-matlab


text:
function v = compute_rbf(x,d,xi, name)

switch name
    case 'abs'
        f = @(x)abs(x);
    case 'gauss'
        sigma = .03;
        f = @(x)exp(-x.^2/(2*sigma^2));
    case 'poly3'
        q = 3;
        f = @(x)abs(x).^q;
    case 'sqrt'
        q = .5;
        f = @(x)abs(x).^q;
    case 'thinplate'
        f = @(x)(x.^2).*log(abs(x)+.001);
end

n = length(xi);
m = length(x);

[Y,X] = meshgrid(x,x);

D = f(X-Y);

% solve for weights
a = pinv(D)*d;

v = f( repmat(xi, [1 m]) - repmat(x', [n 1]) ) .* repmat( a', [n 1] );
v = sum(v, 2);

author: gpeyre

repo: matlab-toolboxes


text:
function [dst,shear_f]=nsst_dec1e(x,shear_parameters,lpfilt)
% This function computes the (local) nonsubsampled shearlet transform as given
% in G. Easley, D. Labate and W. Lim, "Sparse Directional Image Representations
% using the Discrete Shearlet Transform", Appl. Comput. Harmon. Anal. 25 pp.
% 25-46, (2008). This is the more efficient version. Efficiency increases
% if shearing filters (shear_f) are previously stored.
%
% Inputs:
%
% x                       - input image
%
% shear_parameters has the following fields:
%
% shear_parameters.dcomp - a vector such that .dcomp(i) indicates that the
%                          ith decomposition level has 2^decomp(i)
%                          directions. The length of the vector plus 1 is
%                          total the number of decompostions.
%
% shear_parameters.dsize - a vector indicating the local support of the
%                          shearing filter is .dsize(i) for 2^dcomp(i)
%                          directions. This vector is same size as .dcomp.
%
% lpfilt                 - lpfilt is the filter to be used for the Laplacian
%                          Pyramid/ATrous decomposition using the codes
%                          written by Arthur L. Cunha
%
% Output:
%
% dst                    - the cell array containing the discrete shearlet
%                          tranform coefficients
%
% Code contributors: Glenn R. Easley, Demetrio Labate, and Wang-Q Lim.
% Copyright 2011 by Glenn R. Easley. All Rights Reserved.
%

[L,L]=size(x);
level=length(shear_parameters.dcomp);

% LP decomposition
y = atrousdec(x,lpfilt,level);

dst = cell(1,level+1);
dst{1}=y{1}; % assign low-pass coefficients to first decomposition index

shear_f=cell(1,level); % declare cell array containing shearing filters
for i=1:level,
    shear_f{i}=shearing_filters_Myer(shear_parameters.dsize(i),shear_parameters.dcomp(i)).*sqrt(shear_parameters.dsize(i));
    for k=1:2^shear_parameters.dcomp(i),
        dst{i+1}(:,:,k)=conv2(y{i+1},shear_f{i}(:,:,k),'same');
    end
end

author: thomaskuestner


text:
%{
    Demonstration of dimensionality reduction using KPCA.
%}
clc
clear all
close all
addpath(genpath(pwd))

load('.\data\helix.mat', 'data')
kernel = Kernel('type', 'gaussian', 'gamma', 2);
parameter = struct('numComponents', 2, ...
                   'kernelFunc', kernel);
% build a KPCA object
kpca = KernelPCA(parameter);
% train KPCA model
kpca.train(data);

%　mapping data
mappingData = kpca.score;

% Visualization
kplot = KernelPCAVisualization();
% visulize the mapping data
kplot.score(kpca)

author: iqiukp

repo: KPCA-MATLAB

save_path: github-repos/MATLAB/iqiukp-KPCA-MATLAB


text:
% CVX: Keywords for declaring variables and objectives
%   dual      - Declares one or more dual variables.
%   epigraph  - Declares an epigraph variable.
%   hypograph - Declares a hypograph variable.
%   maximize  - Specifies a concave (or affine) objective to be maximized.
%   minimize  - Specifies a convex (or affine) objective to be maximized.
%   subject   - Implements the "subject to" keyword.
%   variable  - Declares a single CVX variable with optional matrix structure.
%   variables - Declares one or more CVX variables.

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST