text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_funcNoise(Temp,Event,handles)
handles=guidata(handles.MU_matrix_display);

[Type,ok] = listdlg('ListString',{'Gaussian white noise','Salt & Pepper','Multiplicative noise (speckle)'}, ...
                         'SelectionMode','single',...
                         'PromptString','Noise Type',...
                         'Name','Type');
if ok==0
    warndlg('Adding noise is cancelled.');
    return;
end

switch Type
    case 1 % gaussian
        f = inputdlg({'Mean:', 'Variance:'},'Gaussian White Noise',1,{'0' ,'1'});
        if isempty(f)
            warndlg('Adding noise is cancelled.');
            return;
        end
    case 2 % salt & pepper
        f = inputdlg({'Noise percent (%):'},'Salt & Pepper',1,{'5'});
        if isempty(f)
            warndlg('Adding noise is cancelled.');
            return;
        end
    case 3 % speckle
        f = inputdlg({'Noise variance:'},'Multiplicative Noise',1,{'0.04'});
        if isempty(f)
            warndlg('Adding noise is cancelled.');
            return;
        end
end

choice = questdlg('Apply to all slices?','All Slices','No','Yes','No');
if isempty(choice)
    warndlg('Adding noise is cancelled.');
    return;
end
% Handle response
switch choice
    case 'No'
        TMatrix = handles.BMatrix;
    case 'Yes'
        TMatrix = handles.TMatrix;
end

TMatrix = double(TMatrix); % convert to double

switch Type
    case 1
        TMatrix = TMatrix + str2double(f{1}) + sqrt(str2double(f{2})) * randn(size(TMatrix));
    case 2
        [row,col,layer] = size(TMatrix);
        maxT=max(TMatrix(:));
        minT=min(TMatrix(:));
        numPx = str2double(f{1}) * row * col *layer /100;
        Idx = randi(row * col * layer, [round(numPx), 1]);
        TMatrix(Idx(1:floor(end/2))) = maxT;
        TMatrix(Idx(ceil(end/2):end)) = minT;
    case 3
        TMatrix = TMatrix + ((2*sqrt(3)*sqrt(str2double(f{1})))*rand(size(TMatrix)) + (-sqrt(3)*sqrt(str2double(f{1})))) .* TMatrix;
end

TMatrix = cast(TMatrix, class(handles.BMatrix)); % convert back to whatever original

switch choice
    case 'No'
        handles.TMatrix(:,:,handles.V.Slice) = TMatrix;
    case 'Yes'
        handles.TMatrix = TMatrix;
end

MergeM=get(handles.Matrix_name_edit,'String');
set(handles.Matrix_name_edit,'String',[MergeM '_noi']);

% update current display matrix
handles=MU_update_image(handles.Matrix_display_axes,{handles.TMatrix,handles.Mask},handles,0);
guidata(handles.MU_matrix_display, handles);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function ori = project(ori,varargin)
% projects orientation to a fundamental region
%
% Syntax
%   ori = project(ori,f)
%
% Input
%  ori     - @orientation
%  ori_ref - reference @rotation
%  f       - @fibre
%
% Output
%  phi1, Phi, phi2 - Euler angles
%

if isa(varargin{1},'fibre') % project to fibre

  %f = varargin{1};

  %h = f.h; r = f.r;
  h = varargin{1}; r = varargin{2};
  rotSym = ori.CS.properGroup.rot;

  d = dot_outer(times(rotSym,h,1), times(inv(ori),r,0),'noSymmetry').';
  [~,idSym] = max(d,[],2);
  ori = times(ori, rotSym(idSym),0);

else

end

author: mtex-toolbox

repo: mtex

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


% =========================================================================
% This code is part of the Matlab-based toolbox
% LagLDDDM - A Lagrangian Gauss--Newton--Krylov Solver for Mass- and
%                        Intensity-Preserving Diffeomorphic Image Registration
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM
%
% Setup NIREP problem (data needs to be obtained
% separately)
%
% =========================================================================

dataset='na03';
example = ['3D-nirep-',dataset];
checkSetupDataFile; if OK, return; end;

% set view options and interpolation options
[viewer,viewOptn] = viewImage('reset','viewImage','imgmontage','colormap','gray(256)','direction','-zyx');

% setup interpolation scheme
imgOptn  = {'imgModel','splineInterMex','regularizer','moments','theta',.01};

% setup transformation used in the parametric part
traOptn  = {'trafo','affine3D'};

% setup distance measure
disOptn  = {'distance','SSD'};

% initialize the regularizer for the non-parametric part
regOptn = {'regularizer','mbHyperElastic','alpha',1,'alphaLength',1,'alphaArea',.1,'alphaVolume',2};

FAIRmessage(mfilename)
load('na01-128x150x128');
load('na01-128x150x128-labels');
load([dataset,'-128x150x128']);
load([dataset,'-128x150x128-labels']);
dataT = double(dataT);
dataR = double(dataR);
dataTl = double(dataTl);
dataRl = double(dataRl);
m = size(dataT);

dataT = 256.*dataT;
dataR = 256.*dataR;
max(dataT(:))

%omega = [0,m(1),0,m(2),0,m(3)];
omega = [0,20,0,23.4375,0,20];
ML = getMultilevel({dataT,dataR},omega,m,'fig',2);
save(outfile,'dataT','dataR','dataTl','dataRl','omega','m','ML');
save(outfile,'-append','viewOptn','imgOptn','traOptn','disOptn','regOptn');
checkSetupDataFile;

% xc       = getCellCenteredGrid(omega,m);
% viewData = @(I) viewImage(imgModel(I,omega,xc),omega,m);

% FAIRfigure(1,'figname',mfilename); clf;
% subplot(1,2,1); viewData(dataT); title('template');
% subplot(1,2,2); viewData(dataR); title('reference');


text:
close all; clear all;

%% Setting of the problem
global s
pde = fracLapdata6;
option.theta = 0.3;
option.estType = 'star';
option.maxIt = 18;
option.maxN = 5e4;
option.solver = 'mg';
option.tol = 1e-6;
[node,elem] = squaremesh([-1,1,-1,1],0.5);
bdFlag = setboundary(node,elem,'Dirichlet');

% %% s = 0.2
% s = 0.2;
% afemfracLap(node,elem,pde,bdFlag,option);
%
% %% s = 0.4
% s = 0.4;
% afemfracLap(node,elem,pde,bdFlag,option);

%% s = 0.6
s = 0.6;
afemfracLap(node,elem,pde,bdFlag,option);

%% s = 0.8
s = 0.8;
afemfracLap(node,elem,pde,bdFlag,option);

repo: ifem


text:
    sched = findResource('scheduler','configuration','local');
    pjob  = createParallelJob('configuration', 'local',...
                    'PathDependencies',{'/Users/ecliff/Desktop/lyrics/'});

    % Set the number of workers
    set(pjob, 'MinimumNumberOfWorkers', 3);
    set(pjob, 'MaximumNumberOfWorkers', 3);
    hndl=@systolicLyricsFilter;
    task= createTask(pjob, hndl, 1, {'mjam.wav'});

    submit(pjob);
    waitForState(pjob, 'finished');
    disp('Parallel Job is finished')
    results = getAllOutputArguments(pjob);
%
%  Clean up the job
%
  destroy(pjob);

    [tune, fs] = wavread('mjam.wav');
    disp('Play the original version')
    sound(tune, fs)

    disp('Play the filtered version')
    sound(results{3}, fs)

author: johannesgerer

repo: jburkardt-m


text:
function [ph_out,ph_out_low]=wrap_filt(ph,n_win,alpha,n_pad,low_flag)
%WRAP_FILT Goldstein adaptive and lowpass filtering
%   [ph_out]=wrap_filt(ph,n_win,alpha,n_pad)
%
%   Andy Hooper, June 2006
%
%   Change Log:
%   07/2006 AH: Added zero padding
%   02/2012 AH: Set magnitude to original values

if nargin<4 | isempty(n_pad)
    n_pad=round(n_win*0.25);
end

if nargin<5
    low_flag='n';
end

[n_i,n_j]=size(ph);
n_inc=floor(n_win/2);
n_win_i=ceil(n_i/n_inc)-1;
n_win_j=ceil(n_j/n_inc)-1;

ph_out=zeros(size(ph));
if strcmpi(low_flag,'y')
    ph_out_low=ph_out;
else
    ph_out_low=[];
end
x=[1:n_win/2];
[X,Y]=meshgrid(x,x);
X=X+Y;
wind_func=[X,fliplr(X)];
wind_func=[wind_func;flipud(wind_func)];

ph(isnan(ph))=0;
B=gausswin(7)*gausswin(7)';
ph_bit=zeros(n_win+n_pad);

L=ifftshift(gausswin(n_win+n_pad,16)*gausswin(n_win+n_pad,16)');

for ix1=1:n_win_i
    wf=wind_func;
    i1=(ix1-1)*n_inc+1;
    i2=i1+n_win-1;
    if i2>n_i
        i_shift=i2-n_i;
        i2=n_i;
        i1=n_i-n_win+1;
        wf=[zeros(i_shift,n_win);wf(1:n_win-i_shift,:)];
    end
    for ix2=1:n_win_j
        wf2=wf;
        j1=(ix2-1)*n_inc+1;
        j2=j1+n_win-1;
        if j2>n_j
           j_shift=j2-n_j;
           j2=n_j;
           j1=n_j-n_win+1;
           wf2=[zeros(n_win,j_shift),wf2(:,1:n_win-j_shift)];
        end
        ph_bit(1:n_win,1:n_win)=ph(i1:i2,j1:j2);
        ph_fft=fft2(ph_bit);
        H=abs(ph_fft);
        H=ifftshift(filter2(B,fftshift(H))); % smooth response
        meanH=median(H(:));
        if meanH~=0
            H=H/meanH;
        end
        H=H.^alpha;
        ph_filt=ifft2(ph_fft.*H);
        ph_filt=ph_filt(1:n_win,1:n_win).*wf2;
        if strcmpi(low_flag,'y')
            ph_filt_low=ifft2(ph_fft.*L);
            ph_filt_low=ph_filt_low(1:n_win,1:n_win).*wf2;
        end
        if isnan(ph_filt(1,1))
            disp('filtered phase contains NaNs in goldstein_filt')
            keyboard
        end
        ph_out(i1:i2,j1:j2)=ph_out(i1:i2,j1:j2)+ph_filt;
        if strcmpi(low_flag,'y')
            ph_out_low(i1:i2,j1:j2)=ph_out_low(i1:i2,j1:j2)+ph_filt_low;
        end
    end
end

ph_out=abs(ph).*exp(j*angle(ph_out)); % reset magnitude
if strcmpi(low_flag,'y')
    ph_out_low=abs(ph).*exp(j*angle(ph_out_low)); % reset magnitude
end

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
% STD_FILEINFO - Check uniform channel distribution across datasets
%
% Usage:
%   >> [struct filepresent] = std_fileinfo(ALLEEG);
% Inputs:
%   ALLEEG    - EEGLAB ALLEEG structure
%
% Outputs:
%   struct      - structure of the same length as the ALLEEG variable
%                 containing all the fields in the datafiles
%   filepresent - array of 0 and 1 indicating if the file is present for
%                 each dataset
%
% Authors: Arnaud Delorme, SCCN/UCSD, CERCO/CNRS, 2010-

% Copyright (C) Arnaud Delorme
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [tmpstructout compinds filepresent] = std_fileinfo( ALLEEG, filetype );

    firstpass   = 1;
    notequal = 0;
    compinds = {};
    tmpstructout = [];
    filepresent = zeros(1,length(ALLEEG));
    for dat = 1:length(ALLEEG)

        filename = fullfile(ALLEEG(dat).filepath, [ ALLEEG(dat).filename(1:end-3) filetype ]);
        thisfilepresent = exist(filename);
        if thisfilepresent && firstpass == 1
            %fprintf('Files of type "%s" detected, checking...',  filetype);
        elseif firstpass == 1
            notequal = 1;
        end
        firstpass = 0;
        filepresent(dat) = thisfilepresent;

        if filepresent(dat)
            try
                warning('off', 'MATLAB:load:variableNotFound');
                tmptmpstruct = load( '-mat', filename, 'times', 'freqs', 'parameters', 'labels', 'chanlabels' );
                warning('on', 'MATLAB:load:variableNotFound');
            catch
                passall = 0;
                fprintf(' Error\n');
                break;
            end

            % rename chanlabels and store structure
            if isfield(tmptmpstruct, 'chanlabels')
                tmptmpstruct.labels = tmptmpstruct.chanlabels;
                tmptmpstruct = rmfield(tmptmpstruct, 'chanlabels');
            end

            if filetype(1) ~= 'd' % ICA components
                allvars = whos('-file', filename);
                tmpinds = [];
                for cind = 1:length(allvars)
                    str = allvars(cind).name(5:end);
                    ind_ = find(str == '_');
                    if ~isempty(ind_), str(ind_:end) = []; end
                    tmpinds = [ tmpinds str2num(str) ];
                end
                compinds(dat) = { unique(tmpinds) };
            elseif ~isfield(tmptmpstruct, 'labels')
                allvars = whos('-file', filename);
                allvarnames = { allvars.name };
                tmptmpstruct.labels = allvarnames(strmatch('chan', allvarnames));
            end
            try,
                tmpstruct(dat) = tmptmpstruct;
            catch,
                passall = 0;
                break;
            end
        end
    end
    if exist('tmpstruct') == 1
        tmpstructout = tmpstruct;
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function h = som_barplane(varargin)

%SOM_BARPLANE Visualize the map prototype vectors as bar charts
%
% h = som_barplane(lattice, msize, data, [color], [scaling], [gap], [pos])
% h = som_barplane(topol, data, [color], [scaling], [gap], [pos])
%
%  som_barplane('hexa',[5 5], rand(25,4), jet(4))
%  som_barplane(sM, sM.codebook,'none')
%
%  Input and output argumetns ([]'s are optional):
%   lattice   (string) grid 'hexa' or 'rect'
%   msize     (vector) size 1x2, defines the map grid size msize, M=prod(msize)
%             (matrix) size Mx2, gives explicit coordinates for each node:
%                      in this case the first argument does not matter.
%   topol     (struct) map or topology struct
%   data      (matrix) size Mxd, each row defines heights of the bars
%   [color]   (matrix) size dx3, of RGB triples. The rows define colors
%                      for each bar in a node. Default is hsv(d). A ColorSpec or
%             (string) A ColorSpec or 'none' gives each bar the same color.
%   [scaling] (string) 'none', 'unitwise' or 'varwise'. The scaling
%                      mode for the values. Default is 'varwise'.
%   [gap]     (scalar) Defines the gap between bars, limits: 0 <= gap <= 1
%                      where 0=no gap, 1=bars are thin lines. Default is 0.25.
%   [pos]     (vector) 1x2 vector defines the position of origin.
%                      Default is [1 1].
%
%   h         (scalar) the object handle to the PATCH object
%
% Axis are set as in SOM_CPLANE.
%
% For more help, try 'type som_barplane' or check out online documentation.
% See also SOM_CPLANE, SOM_PLOTPLANE, SOM_PIEPLANE.

%%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% som_barplane
%
% PURPOSE
%
% Visualizes the map prototype vectors as bar charts.
%
% SYNTAX
%
%  h = som_barplane(topol, data)
%  h = som_barplane(lattice, msize, data)
%  h = som_barplane(..., color)
%  h = som_barplane(..., color, scaling)
%  h = som_barplane(..., color, scaling, gap)
%  h = som_barplane(..., color, scaling, gap, pos)
%
% DESCRIPTION
%
% Visualizes the map prototype vectors as bar charts.
%
% REQUIRED INPUT ARGUMENTS
%
% lattice     The basic shape of the map units
%    (string) 'hexa' or 'rect' positions the bar charts according to
%             hexagonal or rectangular map lattice
%
% msize       The size of the map grid
%    (vector) [n1 n2] vector defines the map size (height: n1 units widht: n2
%             units, total: M=n1xn2 units). The units will be placed to their
%             topological locations in order to form a uniform hexagonal or
%             rectangular grid.
%    (matrix) Mx2 matrix defines arbitary coordinates for the N units. In
%             this case the argument 'lattice' has no effect
%
% topol    Topology of the map grid
%
%   (struct) map or topology struct from which the topology is taken
%
% data        The data to use when constructing the bar charts.
%             Typically, the map codebook or some of its components.
%    (matrix) Mxd matrix. A row defines heights of the bars.
%
% OPTIONAL INPUT ARGUMENTS
%
% Note: if unspecified or given an empty value ('' or []), default
% values are used for optional input arguments.
%
% color       The color of the bars in each pie
%    (ColorSpec) or (string) 'none' gives the same color for each slice.
%    (matrix) dx3 matrix assigns an RGB color determined by the dth row of
%             the matrix to the dth bar (variable) in each bar plot.
%             Default value is hsv(d).
%
% scaling     How to scale the values
%    (string) 'none', 'unitwise' or 'varwise'. This determines the
%             scaling of codebook values when drawing the bars.
%
%             'none' don't scale at all. The bars are not limited
%              to remain inside he units' area: That is, if value of
%              some variable exceeds [-.625,.625] for 'rect' (and
%              in "worst case" [-.5,-.5] for 'hexa') the bars may
%              overlap other units.
%
%              Base line (zero value line)
%              - is in the middle of the unit if data (codebook) contains both
%                negative and positive values (or is completely zero).
%              - is in the top the unit if data (codebook) contains only
%                non-positive values (everything <=0).
%              - is in the bottom the unit if data (codebook) contains only
%                non-negative values (everything >=0).
%
%             'varwise' scales values so that each variable is scaled separately
%              so that when it gets its overall maximum value, the
%              corresponding bar gets maximum range and for minimum value
%              it gets the minimum range. Baseline: see scaling 'none'
%              This is the default.
%
%             'unitwise' scales values in each unit individually so that the
%              bars for variables having minimum and maximum values have minimum
%              and maximum range inside each unit, respectively.
%              In this case the zero value line may move depending on the values.
%

% gap         The gap between bars
%    (scalar) 0: no gap: bars are glued together
%             ... default value is 0.25
%             1: maximum gap: bars are thin lines
%
% pos         Position of origin
%    (vector) size 1x2. This is meant for drawing the plane in arbitrary
%             location in a figure. Note the operation: if this argument is
%             given, the axis limits setting part in the routine is skipped and
%             the limits setting will be left to be done by MATLAB's defaults.
%             Default is [1 1].
%
% OUTPUT ARGUMENTS
%
%  h (scalar) handle to the created patch object
%
% OBJECT TAGS
%
%  One object handle is returned: field Tag is set to 'planeBar'
%
% FEATURES
%
%  - The colors are fixed: changing colormap in the figure (see help
%    colormap) will not change the coloring of the bars.
%
% EXAMPLES
%
% %%% Create the data and make a map
%
% data=rand(100,5); map=som_make(data);
%
% %%% Create a 'jet' colormap that has as many rows as the data has variables
%
% colors=jet(5);
%
% %%% Draw bars
%
% som_barplane(map.topol.lattice, map.topol.msize, map.codebook, colors);
% or som_barplane(map.topol, map.codebook, colors);
% or som_barplane(map, map.codebook, colors);
%
% %%% Draw the bars so that the gap between the bars is bigger and all
%     bars are black
%
% som_barplane(map, map.codebook, 'k', '', 0.6);
%
% SEE ALSO
%
% som_cplane     Visualize a 2D component plane, u-matrix or color plane
% som_plotplane  Visualize the map prototype vectors as line graphs
% som_pieplane   Visualize the map prototype vectors as pie charts

% Copyright (c) 1999-2000 by the SOM toolbox programming team.
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 2.0beta Juha P 110599, Johan 140799, juuso 151199 140300 070600

%%% Check & Init arguments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nargin, lattice, msize, data, color, scaling, gap, pos] = vis_planeGetArgs(varargin{:});
error(nargchk(3, 7, nargin))   % check that no. of input args is correct

% Check pos

if nargin < 7 || isempty(pos)
  pos=NaN;                            % default value for pos (no translation)
elseif ~vis_valuetype(pos,{'1x2'})
  error('Position of origin has to be given as an 1x2 vector');
end

% Check gap

if nargin < 6 || isempty(gap),
  gap=0.25;                           % default value for gap
elseif ~vis_valuetype(gap, {'1x1'}),
  error('Gap value must be scalar.');
elseif ~(gap >= 0 && gap<=1)
  error('Gap value must be in interval [0,1].')
end

% Check scaling

if nargin < 5 || isempty(scaling),
  scaling='varwise';
elseif ~vis_valuetype(scaling,{'string'}) || ...
      ~any(strcmp(scaling,{'none','unitwise','varwise'})),
  error('scaling sholud be ''none'', ''unitwise'' or ''varwise''.');
end

% Check msize

if ~vis_valuetype(msize,{'1x2','nx2'}),
  error('msize has to be 1x2 grid size vector or a Nx2 coordinate matrix.');
end

% Check data

if ~isnumeric(data),
  error('Data matrix has to be numeric.');
elseif length(size((data)))>2
  error('Data matrix has too many dimensions!');
else
  d=size(data,2);
  N=size(data,1);
end

s=.8;                                 % patch size scaling factor

switch scaling,
case 'none'
  % no scaling: don't scale
  % Check data max and min values
  positive=any(data(:)>0); negative=any(data(:)<0);
  if (positive && negative) || (~positive && ~negative),
    % Data contains both negative and positive values (or is
    % completely zero) baseline to centre
    zeroline='zero';
  elseif positive && ~negative
    % Data contains only positive values: baseline to bottom
    zeroline='bottom';
  elseif ~positive && negative
    % Data contains only negative values: baseline to top
    zeroline='top';
  end
case 'unitwise'
  % scale the variables so that the bar for variable with the maximum
  % value in the unit spans to the upper edge of the unit
  % and the bar for the variable with minimum value spans to the lower edge,
  % respectively.
  zeroline='moving';
 case 'varwise'
  % Check data max and min values
  positive=any(data(:)>0); negative=any(data(:)<0);
  if (positive && negative) || (~positive && ~negative),
    % Data contains both negative and positive values (or is
    % completely zero) baseline to
    % centre, scale data so that it doesn't overflow
    data=data./repmat(max(abs([max(data); min(data)])),N,1)*.5;
    zeroline='zero';
  elseif positive && ~negative
    % Data contains only positive values: baseline to
    % bottom, scale data so that it doesn't overflow
    data=data./repmat(max(abs([max(data); min(data)])),N,1)*.5;
    zeroline='bottom';
  elseif ~positive && negative
    % Data contains only negative values: baseline to
    % top, scale data so that it doesn't overflow
    zeroline='top';
    data=data./repmat(max(abs([max(data); min(data)])),N,1)*.5;
  end
otherwise
  error('Unknown scaling mode?');
end

for i=1:N,                            % calculate patch coordinates for
  v=data(i,:);
  [nx,ny]=vis_barpatch(v,gap,zeroline); % bars
  barx(:,(1+(i-1)*d):(i*d))=s*nx;
  bary(:,(1+(i-1)*d):(i*d))=s*ny;
end
l=size(barx,1);

if size(msize,1) == 1,
  xdim=msize(2);
  ydim=msize(1);
  if xdim*ydim~=N
    error('Data matrix has wrong size.');
  else
    y=reshape(repmat(1:ydim,d,1),1,d*ydim); y=repmat(repmat(y,l,1),1,xdim);
    x=reshape(repmat(1:xdim,l*ydim*d,1),l,N*d);
  end
else
  x=reshape(repmat(msize(:,1),1,l*d)',l,d*N);
  y=reshape(repmat(msize(:,2),1,l*d)',l,d*N);
  if N ~= size(msize,1),
    error('Data matrix has wrong size.');
  else
    lattice='rect';
    if isnan(pos),
      pos=[0 0];
    end
  end
end

% Check lattice

if ~ischar(lattice)
  error('Invalid lattice.');
end

switch lattice
case {'hexa','rect'}
  pos=pos-1;
otherwise
  error([ 'Lattice' lattice ' not implemented!']);
end

% Check color
% C_FLAG is for color 'none'

if nargin < 4 || isempty(color)
  color=hsv(d);                       % default n hsv colors
end
if ~vis_valuetype(color, {[d 3],'nx3rgb'},'all') && ...
  ~vis_valuetype(color,{'colorstyle','1x3rgb'})
error('The color matrix has wrong size or has invalid values.');
elseif ischar(color) && strcmp(color,'none')
  C_FLAG=1;
  color='w';
else
  C_FLAG=0;
end

%% Action %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Making lattice.
% Command view([0 90]) shows the map in 2D properly oriented

switch lattice
  case 'hexa'
   t=find(rem(y(1,:),2));                        % move even rows by .5
   x(:,t)=x(:,t)-.5;
   x=x+barx+.5;
   y=y+bary;
  case 'rect'
   x=x+barx;
   y=y+bary;
end

% NB: The coordinates in hexa are not uniform in order to get even
% y-coordinates for the nodes. This is handled by setting _axis scaling_
% so that the hexa-nodes look like uniform hexagonals. See
% vis_PlaneAxisProperties

if ~isnan(pos)
  x=x+pos(1);y=y+pos(2);               % move upper left corner
end                                    % to pos

%% Set axes properties

ax=newplot;                            % get current axis
vis_PlaneAxisProperties(ax,lattice, msize, pos);

%% Rearrange dx3 color matrix

if ~ischar(color) && size(color,1)~=1,
  color=reshape(repmat(color,N,1),[1 N*d 3]);
end

%% Draw the plane!

if isnumeric(color),
  % explicit color settings by RGB-triplets won't work with
  % patch in 'painters' mode, unless there only a single triplet
  si = size(color);
  if length(si)~=2 || any(si==[1 3]), set(gcf,'renderer','zbuffer'); end
end

h_=patch(x,y,color);

if C_FLAG
  set(h_,'FaceColor','none');
end

set(h_,'Tag','planeBar');              % tag the object

%%% Build output %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargout>0, h=h_; end                % Set h only if
                                       % there really is output

%%% Subfunctions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xcoord,ycoord]=vis_barpatch(y,gap,zeroline)

x = length(y);
d = gap/(2*(x-1)+2);
step= -.5:1/x:.5;

miny=min(y);
maxy=max(y);

switch(zeroline)
 case 'moving'
  if miny < 0
    if maxy > 0
      zl = .5 - (abs(miny)/(maxy-miny));    %reverse mode
      y= .5 - ((y-miny*ones(1,x))./(maxy-miny));
    else
      zl = -.5;
      y=-.5+abs(y./miny);
    end
  else
    zl = .5;    %reverse mode
    y=.5-y./maxy;
  end
 case 'moveNotScale'
  if miny < 0
    if maxy > 0
      zl = 0.5+miny;
      y = zl - y;
    else
      zl=-.5;
      y=-.5+abs(y);
    end
  else
    zl=.5;
    y =.5-y;
  end
 case 'zero'
  zl=0; y=zl-y;
 case 'top'
  zl=-.5; y=zl-2*y;
 case 'bottom'
  zl=.5; y=zl-2*y;
end

for i=1:x
  xcoord(:,i) = [d+step(i);d+step(i);step(i+1)-d;step(i+1)-d;d+step(i)];
  ycoord(:,i) = [zl;y(i);y(i);zl;zl];
end

author: ilarinieminen

repo: SOM-Toolbox


text:
function pde = HodgeLaplacian3Edata1
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

pde.f = @f;
pde.exactu = @exactu;
pde.exactsigma = @exactsigma;
pde.gu = @exactu;
pde.gsigma = @exactsigma;
pde.gun = @gun;
pde.curlu = [0 0 0];
pde.gcurlu = 0;

    function s = f(p)
    x = p(:,1); y = p(:,2);
    s = [cos(x), -sin(y), zeros(size(p,1),1)];
    end

    function s = exactu(p)
    x = p(:,1); y = p(:,2); z = p(:,3);
    s = [cos(x), -sin(y), 2*z+1];
    end

    function s = exactsigma(p)
    x = p(:,1); y = p(:,2);
    s = sin(x)+cos(y)-2;
    end

    function f = gun(p) % for unit cube [0,1]^3
    f = zeros(size(p,1),1);
    x = p(:,1); y = p(:,2); z = p(:,3);
    u = exactu(p);
    leftbd = (abs(x)<eps);    % n = (-1,0,0)
    f(leftbd) = - u(leftbd,1);
    rightbd = (abs(x-1)<eps); % n = (1,0,0)
    f(rightbd) = u(rightbd,1);
    frontbd = (abs(y)<eps);   % n = (0,-1,0)
    f(frontbd) = -u(frontbd,2);
    backbd = (abs(y-1)<eps);  % n = (0,1,0)
    f(backbd) = u(backbd,2);
    topbd = (abs(z-1)<eps);   % n = (0,0,1)
    f(topbd) = u(topbd,3);
    bottombd = (abs(z)<eps);  % n = (0,0,-1)
    f(bottombd) = - u(bottombd,3);
    end

end

repo: Aero-matlab


text:
%  plots the curvature of an estimated parameter at the model
%
%  ::
%
%    [h,legend_]=curvature(xname,db)
%
%  Args:
%
%     - **xname** [char]: name of the parameter to plot
%     - **db** [struct]: structure containing the various parameters. Each
%       parameter field is itself a structure with the following
%
%       - **tex_name** [char]: name of the parameter as it should appear in the
%         title of the plot
%       - **mode** [scalar]: value of the parameter at the mode
%       - **log_post_mode** [scalar]: value of the log posterior mode
%       - **log_lik_mode** [scalar]: value of the log likelihood at the mode
%       - **x** [vector]: x-axis values
%       - **log_post** [vector]: value of the log-posterior for each value of x
%       - **log_lik** [vector]: value of the log-likelihood for each value of x
%
%  Returns:
%     :
%
%     - **h** [handle]: handle for the plot
%     - **legend_** [cellstr]: names of the lines in the plot
%     - **tex_name** [char]: name of the parameter
%
%

author: jmaih


text:
% demo script for regression with HKL
clear all

% fixing the seed of the random generators
seed=1;
randn('state',seed);
rand('state',seed);

% toy example characteristics
p = 1024;           % total number of variables (used to generate a Wishart distribution)
psub = 32;          % kept number of variables = dimension of the problem
n = 256;            % number of observations
s = 4;              % number of relevant variables
noise_std = .2;		% standard deviation of noise
proptrain = .5;     % proportion of data kept for training (the rest is used for testing)

% generate random covariance matrix from a Wishart distribution
Sigma_sqrt = randn(p,p);
Sigma = Sigma_sqrt' * Sigma_sqrt;

% normalize to unit trace and sample
diagonal = diag(Sigma);
Sigma = diag( 1./diagonal.^.5) * Sigma * diag( 1./diagonal.^.5);
Sigma_sqrt =   Sigma_sqrt * diag( 1./diagonal.^.5);
X = randn(n,p) * Sigma_sqrt;

X = X(:,1:psub);
p=psub;

% generate nonlinear function of X as the sum of all cross-products
J =  1:s;    % select the first s variables
Y = zeros(n,1);
for i=1:s
	for j=1:i-1
		Y = Y + X(:,J(i)) .* X(:,J(j));
	end
end
% normalize to unit standard deviation
Y = Y / std(Y);

% add some noise with known standard deviation
Y =  Y + randn(n,1) * noise_std;

% split data in two groups
ntrain = round(n*proptrain);
ntest = n - ntrain;
rp = randperm(n);
trainset = rp(1:ntrain);
testset  = rp(ntrain+1:end);

% select some regularization parameterst lambdas (from large to small)
lambdas = 10.^[2:-.5:-8];

% HKL with hermite polynomial decomposition
[outputs,model,accuracies,accuracies_kfold,bestlambda_same] = hkl_kfold(10,'same',X(trainset,:),Y(trainset),lambdas,'square','hermite',[ .5 3 .1 4 ],...
	'Xtest',X(testset,:),'Ytest',Y(testset),'maxactive',300,'memory_cache',2e9,'display',0);

estimated_best_error_same = accuracies.testing_error(bestlambda_same)
best_error_same = min(accuracies.testing_error)

[outputs,model,accuracies,accuracies_kfold,bestlambda_scaled] = hkl_kfold(10,'scaled',X(trainset,:),Y(trainset),lambdas,'square','hermite',[ .5 3 .1 4 ],...
	'Xtest',X(testset,:),'Ytest',Y(testset),'maxactive',300,'memory_cache',2e9,'display',0);

estimated_best_error_scaled = accuracies.testing_error(bestlambda_scaled)
estimated_best_error_overfit = min(accuracies.testing_error)

repo: REST


text:
function [U, S] = pca(X)
%PCA Run principal component analysis on the dataset X
%   [U, S, X] = pca(X) computes eigenvectors of the covariance matrix of X
%   Returns the eigenvectors U, the eigenvalues (on diagonal) in S
%

% Useful values
[m, n] = size(X);

% You need to return the following variables correctly.
U = zeros(n);
S = zeros(n);

% ====================== YOUR CODE HERE ======================
% Instructions: You should first compute the covariance matrix. Then, you
%               should use the "svd" function to compute the eigenvectors
%               and eigenvalues of the covariance matrix.
%
% Note: When computing the covariance matrix, remember to divide by m (the
%       number of examples).
%

Sigma = X'*X/m;

[U,S] = svd(Sigma);

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

end

repo: MachineLearningEx


text:
function [mi2] = m22mi2(m2)
% Convert area from square meters to square miles.
% Chad A. Greene 2012
mi2 = m2*3.861021585425E-7;

author: Sable

repo: mcbench-benchmarks


text:
classdef CEC2013_F5 < PROBLEM
% <single> <real> <large>
% 7-nonseparable, 1-separable shifted and rotated Rastrigin's function

%------------------------------- Reference --------------------------------
% X. Li, K. Tang, M. N. Omidvar, Z. Yang, and K. Qin, Benchmark functions
% for the CEC'2013 special session and competition on large-scale global
% optimization, RMIT University, Australia, 2013.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        Xopt;	% Optimal decision vector
        R25;    % Rotation matrices
        R50;
        R100;
        p;      % Rank of decision variables
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2013.mat'),'Data');
            obj.Xopt = Data{5}.xopt;
            obj.R25  = Data{5}.R25;
            obj.R50  = Data{5}.R50;
            obj.R100 = Data{5}.R100;
            obj.p    = Data{5}.p;
            obj.M    = 1;
            obj.D    = 1000;
            obj.lower    = zeros(1,obj.D) - 5;
            obj.upper    = zeros(1,obj.D) + 5;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            S = [50 25 25 100 50 25 25 700];
            W = [1.81e-1 9.08e3 2.43e1 1.86e-6 1.77e4 2.82e-4 1.53e-2 1];
            PopDec = PopDec - repmat(obj.Xopt,size(PopDec,1),1);
            PopObj = zeros(size(PopDec,1),1);
            for i = 1 : length(S)
                loc = obj.p(sum(S(1:i-1))+1:sum(S(1:i)));
                switch S(i)
                    case 25
                        PopDec(:,loc) = PopDec(:,loc)*obj.R25;
                    case 50
                        PopDec(:,loc) = PopDec(:,loc)*obj.R50;
                    case 100
                        PopDec(:,loc) = PopDec(:,loc)*obj.R100;
                end
                PopObj = PopObj + W(i)*Rastrigin(Tdiag(Tasy(Tosz(PopDec(:,loc)),0.2),10));
            end
        end
    end
end

function F = Rastrigin(X)
    F = sum(X.^2-10*cos(2*pi*X)+10,2);
end

function Z = Tosz(X)
    X1 = zeros(size(X));
    X1(X~=0) = log(abs(X(X~=0)));
    C1 = zeros(size(X)) + 5.5;
    C1(X>0) = 10;
    C2 = zeros(size(X)) + 3.1;
    C2(X>0) = 7.9;
    Z = sign(X).*exp(X1+0.049*(sin(C1.*X1)+sin(C2.*X1)));
end

function Z = Tasy(X,beta)
    Z = X.^(1+repmat(beta*linspace(0,1,size(X,2)),size(X,1),1).*sqrt(X));
    Z(X<=0) = X(X<=0);
end

function Z = Tdiag(X,alpha)
    Z = X.*repmat(sqrt(alpha).^linspace(0,1,size(X,2)),size(X,1),1);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%% DEMO_febio_0023_beam_gravity_body_load
% Below is a demonstration for:
%
% * Building geometry for a beam with hexahedral elements
% * Defining the boundary conditions
% * Coding the febio structure
% * Running the model
% * Importing and visualizing the displacement results

%% Keywords
%
% * febio_spec version 3.0
% * febio, FEBio
% * beam, gravity, body load
% * surface pressure boundary condition
% * hexahedral elements, hex8
% * beam, rectangular
% * static, solid
% * hyperelastic, Ogden
% * displacement logfile

%%

clear; close all; clc;

%% Plot settings
fontSize=20;
faceAlpha1=0.8;
markerSize=40;
markerSize2=20;
lineWidth=3;

%% Control parameters

% Path names
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');

% Defining file names
febioFebFileNamePart='tempModel';
febioFebFileName=fullfile(savePath,[febioFebFileNamePart,'.feb']); %FEB file name
febioLogFileName=[febioFebFileNamePart,'.txt']; %FEBio log file name
febioLogFileName_disp=[febioFebFileNamePart,'_disp_out.txt']; %Log file name for exporting displacements

%Specifying dimensions and number of elements
beamWidth=0.1;
sampleWidth=beamWidth; %Width
sampleLength=6*beamWidth; %Thickness
sampleHeight=beamWidth; %Height
pointSpacings=0.025*ones(1,3); %Desired point spacing between nodes
numElementsWidth=round(sampleWidth/pointSpacings(1)); %Number of elemens in dir 1
numElementsThickness=round(sampleLength/pointSpacings(2)); %Number of elemens in dir 2
numElementsHeight=round(sampleHeight/pointSpacings(3)); %Number of elemens in dir 3

%Material parameter set
c1=1000; %Shear-modulus-like parameter
m1=8; %Material parameter setting degree of non-linearity
k_factor=1e2; %Bulk modulus factor
k=c1*k_factor; %Bulk modulus

%Determine load
sampleVolume=sampleWidth.*sampleLength.*sampleHeight;
sampleWeight=(sampleVolume.*1000);
gravityVector=[0 0 9.81.*sampleWeight];

% FEA control settings
numTimeSteps=10; %Number of time steps desired
max_refs=25; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=6; %Optimum number of iterations
max_retries=5; %Maximum number of retires
dtmin=(1/numTimeSteps)/100; %Minimum time step size
dtmax=1/numTimeSteps; %Maximum time step size
min_residual=1e-40;

%% Creating model geometry and mesh
% A box is created with tri-linear hexahedral (hex8) elements using the
% |hexMeshBox| function. The function offers the boundary faces with
% seperate labels for the top, bottom, left, right, front, and back sides.
% As such these can be used to define boundary conditions on the exterior.

% Create a box with hexahedral elements
beamDimensions=[sampleWidth sampleLength sampleHeight]; %Dimensions
beamElementNumbers=[numElementsWidth numElementsThickness numElementsHeight]; %Number of elements
outputStructType=2; %A structure compatible with mesh view
[meshStruct]=hexMeshBox(beamDimensions,beamElementNumbers,outputStructType);

%Access elements, nodes, and faces from the structure
E=meshStruct.elements; %The elements
V=meshStruct.nodes; %The nodes (vertices)
Fb=meshStruct.facesBoundary; %The boundary faces
Cb=meshStruct.boundaryMarker; %The "colors" or labels for the boundary faces
elementMaterialIndices=ones(size(E,1),1); %Element material indices

%%
% Plotting model boundary surfaces and a cut view

hFig=cFigure;

subplot(1,2,1); hold on;
title('Model boundary surfaces and labels','FontSize',fontSize);
gpatch(Fb,V,Cb,'k',faceAlpha1);
colormap(gjet(6)); icolorbar;
axisGeom(gca,fontSize);

hs=subplot(1,2,2); hold on;
title('Cut view of solid mesh','FontSize',fontSize);
optionStruct.hFig=[hFig hs];
meshView(meshStruct,optionStruct);
axisGeom(gca,fontSize);

drawnow;

%% Defining the boundary conditions
% The visualization of the model boundary shows colors for each side of the
% cube. These labels can be used to define boundary conditions.

%Define supported node set
logicFace=Cb==4; %Logic for current face set
Fr=Fb(logicFace,:); %The current face set
bcSupportList=unique(Fr(:)); %Node set part of selected face

%%
% Visualizing boundary conditions. Markers plotted on the semi-transparent
% model denote the nodes in the various boundary condition lists.

hf=cFigure;
title('Boundary conditions','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;

gpatch(Fb,V,'kw','k',0.5);

hl(1)=plotV(V(bcSupportList,:),'k.','MarkerSize',markerSize);

legend(hl,{'BC support'});

axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Defining the FEBio input structure
% See also |febioStructTemplate| and |febioStruct2xml| and the FEBio user
% manual.

%Get a template with default settings
[febio_spec]=febioStructTemplate;

%febio_spec version
febio_spec.ATTR.version='3.0';

%Module section
febio_spec.Module.ATTR.type='solid';

%Control section
febio_spec.Control.analysis='STATIC';
febio_spec.Control.time_steps=numTimeSteps;
febio_spec.Control.step_size=1/numTimeSteps;
febio_spec.Control.solver.max_refs=max_refs;
febio_spec.Control.solver.max_ups=max_ups;
febio_spec.Control.time_stepper.dtmin=dtmin;
febio_spec.Control.time_stepper.dtmax=dtmax;
febio_spec.Control.time_stepper.max_retries=max_retries;
febio_spec.Control.time_stepper.opt_iter=opt_iter;

%Material section
materialName1='Material1';
febio_spec.Material.material{1}.ATTR.name=materialName1;
febio_spec.Material.material{1}.ATTR.type='Ogden';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.c1=c1;
febio_spec.Material.material{1}.m1=m1;
febio_spec.Material.material{1}.c2=c1;
febio_spec.Material.material{1}.m2=-m1;
febio_spec.Material.material{1}.k=k;

%Mesh section
% -> Nodes
febio_spec.Mesh.Nodes{1}.ATTR.name='nodeSet_all'; %The node set name
febio_spec.Mesh.Nodes{1}.node.ATTR.id=(1:size(V,1))'; %The node id's
febio_spec.Mesh.Nodes{1}.node.VAL=V; %The nodel coordinates

% -> Elements
partName1='Part1';
febio_spec.Mesh.Elements{1}.ATTR.name=partName1; %Name of this part
febio_spec.Mesh.Elements{1}.ATTR.type='hex8'; %Element type
febio_spec.Mesh.Elements{1}.elem.ATTR.id=(1:1:size(E,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=E; %The element matrix

% -> NodeSets
nodeSetName1='bcSupportList';
febio_spec.Mesh.NodeSet{1}.ATTR.name=nodeSetName1;
febio_spec.Mesh.NodeSet{1}.node.ATTR.id=bcSupportList(:);

%MeshDomains section
febio_spec.MeshDomains.SolidDomain.ATTR.name=partName1;
febio_spec.MeshDomains.SolidDomain.ATTR.mat=materialName1;

%Boundary condition section
% -> Fix boundary conditions
febio_spec.Boundary.bc{1}.ATTR.type='fix';
febio_spec.Boundary.bc{1}.ATTR.node_set=nodeSetName1;
febio_spec.Boundary.bc{1}.dofs='x,y,z';

%Loads section
% -> Body load
febio_spec.Loads.body_load{1}.ATTR.type='const';
febio_spec.Loads.body_load{1}.x.ATTR.lc=1;
febio_spec.Loads.body_load{1}.x.VAL=gravityVector(1);
febio_spec.Loads.body_load{1}.y.ATTR.lc=1;
febio_spec.Loads.body_load{1}.y.VAL=gravityVector(2);
febio_spec.Loads.body_load{1}.z.ATTR.lc=1;
febio_spec.Loads.body_load{1}.z.VAL=gravityVector(3);

%LoadData section
% -> load_controller
febio_spec.LoadData.load_controller{1}.ATTR.id=1;
febio_spec.LoadData.load_controller{1}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{1}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{1}.points.point.VAL=[0 0; 1 1];

%Output section
% -> log file
febio_spec.Output.logfile.ATTR.file=febioLogFileName;
febio_spec.Output.logfile.node_data{1}.ATTR.file=febioLogFileName_disp;
febio_spec.Output.logfile.node_data{1}.ATTR.data='ux;uy;uz';
febio_spec.Output.logfile.node_data{1}.ATTR.delim=',';

%% Quick viewing of the FEBio input file structure
% The |febView| function can be used to view the xml structure in a MATLAB
% figure window.

%%
% |febView(febio_spec); %Viewing the febio file|

%% Exporting the FEBio input file
% Exporting the febio_spec structure to an FEBio input file is done using
% the |febioStruct2xml| function.

febioStruct2xml(febio_spec,febioFebFileName); %Exporting to file and domNode

%% Running the FEBio analysis
% To run the analysis defined by the created FEBio input file the
% |runMonitorFEBio| function is used. The input for this function is a
% structure defining job settings e.g. the FEBio input file name. The
% optional output runFlag informs the user if the analysis was run
% succesfully.

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.runMode='external';%'internal';

[runFlag]=runMonitorFEBio(febioAnalysis);%START FEBio NOW!!!!!!!!

%% Import FEBio results

if runFlag==1 %i.e. a succesful run

    %%
    % Importing nodal displacements from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_disp),1,1);

    %Access data
    N_disp_mat=dataStruct.data; %Displacement
    timeVec=dataStruct.time; %Time

    %Create deformed coordinate set
    V_DEF=N_disp_mat+repmat(V,[1 1 size(N_disp_mat,3)]);

    %%
    % Plotting the simulated results using |anim8| to visualize and animate
    % deformations

    DN_magnitude=sqrt(sum(N_disp_mat(:,:,end).^2,2)); %Current displacement magnitude

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    title('Displacement magnitude [mm]','Interpreter','Latex')
    hp=gpatch(Fb,V_DEF(:,:,end),DN_magnitude,'k',1); %Add graphics object to animate
    hp.Marker='.';
    hp.MarkerSize=markerSize2;
    hp.FaceColor='interp';
    gpatch(Fb,V,0.5*ones(1,3),'k',0.25); %A static graphics object

    axisGeom(gca,fontSize);
    colormap(gjet(250)); colorbar;
    caxis([0 max(DN_magnitude)]);
    axis(axisLim(V_DEF)); %Set axis limits statically
    camlight headlight;

    % Set up animation features
    animStruct.Time=timeVec; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments
        DN_magnitude=sqrt(sum(N_disp_mat(:,:,qt).^2,2)); %Current displacement magnitude

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp hp]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData'}; %Properties of objects to animate
        animStruct.Set{qt}={V_DEF(:,:,qt),DN_magnitude}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature
    drawnow;

end

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% Copyright (C) 1999 Paul Kienzle <pkienzle@users.sf.net>
% Copyright (C) 2006 Peter V. Lanspeary, <peter.lanspeary@.adelaide.edu.au>
%
% 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/>.

% -*- texinfo -*-
% @deftypefn  {Function File} {[@var{a}, @var{v}, @var{ref}] =} levinson (@var{acf})
% @deftypefnx {Function File} {[@dots{}] =} levinson (@var{acf}, @var{p})
%
% Use the Durbin-Levinson algorithm to solve:
%    toeplitz(acf(1:p)) * x = -acf(2:p+1).
% The solution [1, x'] is the denominator of an all pole filter
% approximation to the signal x which generated the autocorrelation
% function acf.
%
% acf is the autocorrelation function for lags 0 to p.
% p defaults to length(acf)-1.
% Returns
%   a=[1, x'] the denominator filter coefficients.
%   v= variance of the white noise = square of the numerator constant
%   ref = reflection coefficients = coefficients of the lattice
%         implementation of the filter
% Use freqz(sqrt(v),a) to plot the power spectrum.
%
% REFERENCE
% [1] Steven M. Kay and Stanley Lawrence Marple Jr.:
%   "Spectrum analysis -- a modern perspective",
%   Proceedings of the IEEE, Vol 69, pp 1380-1419, Nov., 1981
% @end deftypefn

% Based on:
%    yulewalker.m
%    Copyright (C) 1995 Friedrich Leisch <Friedrich.Leisch@ci.tuwien.ac.at>
%    GPL license

% FIXME: Matlab doesn't return reflection coefficients and
%        errors in addition to the polynomial a.
% FIXME: What is the difference between aryule, levinson,
%        ac2poly, ac2ar, lpc, etc.?

function [a, v, ref] = oc_levinson (acf, p)

if ( nargin<1 )
    print_usage;
elseif( ~isvector(acf) || length(acf)<2 )
    error( 'levinson: arg 1 (acf) must be vector of length >1\n');
elseif ( nargin>1 && ( ~isscalar(p) || fix(p)~=p ) )
    error( 'levinson: arg 2 (p) must be integer >0\n');
else
    if ((nargin == 1)||(p>=length(acf)))
        p = length(acf) - 1;
    end
    if( size(acf,2)>1 )
        acf=acf(:);
    end      % force a column vector

    if nargout < 3 && p < 100
        % direct solution [O(p^3), but no loops so slightly faster for small p]
        %   Kay & Marple Eqn (2.39)
        R = toeplitz(acf(1:p), conj(acf(1:p)));
        a = R \ -acf(2:p+1);
        a = [ 1, a.' ];
        v = real( a*conj(acf(1:p+1)) );
    else
        % durbin-levinson [O(p^2), so significantly faster for large p]
        %   Kay & Marple Eqns (2.42-2.46)
        ref = zeros(p,1);
        g = -acf(2)/acf(1);
        a = g;
        v = real( ( 1 - g*conj(g)) * acf(1) );
        ref(1) = g;
        for t = 2 : p
            g = -(acf(t+1) + a * acf(t:-1:2)) / v;
            a = [ a+g*conj(a(t-1:-1:1)), g ];
            v = v * ( 1 - real(g*conj(g)) ) ;
            ref(t) = g;
        end
        a = [1, a];
    end
end

end

author: brainstorm-tools


text:
% This code is the implementation of "Multi-scale Guided Image and Video Fusion: A Fast and Efficient Approach"
% Cite this article as:
% Bavirisetti, D.P., Xiao, G., Zhao, J. et al. Circuits Syst Signal Process (2019).
%https://doi.org/10.1007/s00034-019-01131-z
%
% The interface is created by the authors of VIFB.

function img = run_MGFF(imgVI, imgIR, visualization)

    % Guided image filter parameters
    r=9;eps=10^3;

    I1 = double(imread(imgIR.img));
    if size(I1,3)==1
        I1 = repmat(I1,[1,1,3]);
    end

    I2 = double(imread(imgVI.img));

    %% apply multi-scale guided image fusion on source images
    tic
    F = fuse_MGF_RGB(I1, I2, r, eps);
    toc

    fuseimage = im2uint8(F);
    img = fuseimage;

    %% display source images and the fused image
    if visualization == 1
        figure, imshow(uint8(I1), []);
        figure, imshow(uint8(I2),[]);
        figure, imshow((F),[]);
    end
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
% function F = det_F_gold(x1,x2,L_COST,SAMPSON_APPROXIMATION,NORMALIZE)
% Determines the F by iteratively minimizing the geometric error
% Algorithm 11.2 in Hartley & Zisserman, Multiple View Geometry in Computer
% Vision
% Inputs:
%           x1                      3xN coordinates of matched points in image 1(homogeneous)
%           x2                      3xN coordinates of matched points in image 2(homogeneous)
%           L_COST                  1x1 (optional) controls penalization scheme
%                                       for the cost function: L_COST 1 leads to L_1 minimization of
%                                       the average geometric cost (the one mentioned in the book)
%           SAMPSON_APPROXIMATION   1x1 (optional) if enabled,
%                                       approximates the geometric mean with the sampson cost
%           NORMALIZE               1x1 (optional)determines if the algorithm
%                                       should use the normalized points
% Outputs:
%           F                       3x3 the fundamental matrix
%
% Author: Omid Aghazadeh, KTH(Royal Institute of Technology), 2010/05/09
function F = det_F_gold(x1,x2,L_COST,SAMPSON_APPROXIMATION,NORMALIZE)
if sum(size(x1)~=size(x2)), error('size of correspondences do not match!'), end
if size(x1,1) ~= 3, error('invalid points'), end
global MAX_FUN_EVAL MAX_ITER TOL_X TOL_FUN;
if nargin<3, L_COST = 1; end;
if nargin<4, SAMPSON_APPROXIMATION = 0; end
if nargin<5, NORMALIZE = 1; end;
if NORMALIZE
    nmat1 = get_normalization_matrix(x1); % isotropic normalization (translation/scaling)
    nmat2 = get_normalization_matrix(x2);
    x1n = nmat1*x1;
    x2n = nmat2*x2;
else
    nmat1 = eye(3); nmat2 = eye(3); x1n = x1; x2n = x2;
end
x1n = x1n./repmat(x1n(3,:),3,1); x2n = x2n./repmat(x2n(3,:),3,1); % normalizing points so their last coordinate is 1
F_0 = det_F_normalized_8point(x1n,x2n);
%% (ii)

[e,e_prime] = get_epipole(F_0);
e_prime_cross = get_x_cross(e_prime);
P2 = [e_prime_cross*F_0 e_prime];

%% (iii) minimze the cost

if ~ SAMPSON_APPROXIMATION
    [P2] = lsqnonlin(@(p2)costGold(x1n,x2n,p2,L_COST),P2,[],[],optimset('Display','off','TolX',TOL_X,'TolFun',TOL_FUN,'MaxFunEval',MAX_FUN_EVAL,'MaxIter',MAX_ITER,'Algorithm',{'levenberg-marquardt' 0.01}));
else
    [P2] = lsqnonlin(@(p2)costSampson(x1n,x2n,p2),P2,[],[],optimset('Display','off','TolX',TOL_X,'TolFun',TOL_FUN,'MaxFunEval',MAX_FUN_EVAL,'MaxIter',MAX_ITER,'Algorithm',{'levenberg-marquardt' 0.01}));
end

Fhat = get_x_cross(P2(:,4))*P2(:,1:3);

F= nmat2' * Fhat* nmat1; % denormalization

end

%% function scost = costGold(x1,x2,P2,L_COST)
% this is the cost function for the Gold Standard algoritghm. The
% triangulation method is the inhomogeneous one(chapter 12 of the book)
function scost = costGold(x1,x2,P2,L_COST)
Xhat = triangulate(x1,x2,P2,1);
xhat1 = Xhat(1:3,:)./repmat(Xhat(3,:),3,1); % the first camera is assumed to be [I|0]
xhat2 = P2 * Xhat;
xhat2 = xhat2./repmat(xhat2(3,:),3,1);
cost = ((x1(:)-xhat1(:)).^2 + (x2(:)-xhat2(:)).^2);
scost = sqrt(sum(cost))^L_COST;
end

%% function scost = costSampson(x1,x2,P2)
% this is the cost function for the sampson approximation. It implements an
% over-parametrization of F, however the minimal solution can be easily
% integrated here.
function scost = costSampson(x1,x2,P2)
F = get_x_cross(P2(:,4))*P2(:,1:3);
Fx1 = F*x1;
Ftx2 = F'*x2;
num = sum(x2 .* Fx1,1).^2;
denum= sum(Fx1(1:2,:).^2,1) + sum(Ftx2(1:2,:).^2,1);
cost = num./denum;
scost = sqrt(sum(cost));
end


text:
function [xt_prf, sw_prf] = BLanalytical(muw, muo)
%BLANALYTICAL Analytical solution of Buckley-Leverett equation
% Written by Ali A. Eftekhari
% See the license file
%muw = 10e-3;
%muo = 10e-3;
u = 1e-3;
phi = 0.2;
sw_in = 1;
sw_end = 0;
krw = @(sw)(sw.^4);
dkrwdsw = @(sw)(4*sw.^3);
kro = @(sw)((1-sw.^2).*(1-sw).^2);
dkrodsw = @(sw)(-2*sw.*(1-sw).^2-2*(1-sw).*(1-sw.^2));
fw = @(sw)((krw(sw)/muw)./(krw(sw)/muw+kro(sw)/muo));
dfwdsw = @(sw)((dkrwdsw(sw)/muw.*(krw(sw)/muw+kro(sw)/muo)- ...
    (dkrwdsw(sw)/muw+dkrodsw(sw)/muo).*krw(sw)/muw)./ ...
    (krw(sw)/muw+kro(sw)/muo).^2);
s = linspace(0,1,100);
figure(1)
subplot(2,2,1);
plot(s, krw(s), s, kro(s));
F = @(sw)(dfwdsw(sw)-(fw(sw)-fw(sw_end))/(sw-sw_end));
sw_shock = fzero(F, [eps,1]);
subplot(2,2,2);
plot(s, fw(s), [sw_end sw_shock], [fw(sw_end) fw(sw_shock)]);
% plot(s, fw(s), [sw_end sw_shock], [fw(sw_end) fw(sw_shock)]);
s1 = linspace(sw_in, sw_shock, 50);
xt_s1 = u/phi*dfwdsw(s1);
xt_s = u/phi*dfwdsw(s);
xt_shock = u/phi*dfwdsw(sw_shock);
subplot(2,2,3);
plot(xt_s, s, '--', ...
    [xt_s1 xt_shock xt_shock max(xt_s)], [s1 sw_shock sw_end sw_end])
sw_prf = [s1 sw_shock sw_end sw_end];
xt_prf = [xt_s1 xt_shock xt_shock max(xt_s)];
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
classdef AnnotationVideoMarker < handle
    properties (Access = public)
        sample;
        label;
        text;
    end

    methods (Access = public)
        function obj = AnnotationVideoMarker(sample,label,text)
            if nargin ~= 0
                obj.sample = sample;
                obj.label = label;
                obj.text = text;
            end
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
%Note: although this example uses wander-azimuth concept, you should not consider this
%as a suitable wander-azimuth implementaion. This purpose of this example is to illustrate
%how large heading errors can be handled. There is an another example in the toolkit
%presenting a suitable wander-azimuth implementation.

clear;
fclose all;
%PATH='/home/yigiter/Desktop/Dropbox/INS/LargeHeading/mfiles/data/';
PATH='C:\Users\yigiter\Desktop\Dropbox\INS\LargeHeading\mfiles\data\';

%%Specify the imu error model and generate the path
[SenErrDef, IniErrDef, ObsErrDef, ini_pv, wander, dt]=sys_path_largeheading(PATH, 1);

%%Open the files
F_IMU=fopen([PATH 'imu.bin'],'rb');
F_TRU=fopen([PATH 'mnav.bin'],'rb');
F_IMUERR=fopen([PATH 'imuerr.bin'],'rb');
F_OBS_POS=fopen([PATH 'obs_pos.bin'],'rb');
F_RES=fopen([PATH 'res.bin'],'wb');
F_ERR=fopen([PATH 'err.bin'],'wb');

%propagation time for the discretized system
prp_per=round(0.33/dt);

%%Sensor error models
[dAimu, dBimu, Cimu, dDimu, sP0imu]=imu_modTI_v000(SenErrDef);  %errors models in discrete time
[Aimu, Qimu, Rimu]=dc2dc_v000(dAimu, dBimu*dBimu', dDimu*dDimu', dt, 1, []);  %models in continuous time
ximu=zeros(6,1); %imu error states

%%%%Initialize the INS
%PV is provided
pos_g=ini_pv(:,1);
vel_n=ini_pv(:,2); %note path generator is defined for body frame velocity
%Use accelerometer data to initialize the attitude
t=0;
in=0;
imu_avg=zeros(6,1);
imu_data=fread(F_IMU,8,'double'); %garbage
true_val = fread(F_TRU, 10, 'double');
imu_err=fread(F_IMUERR, 7, 'double');
avg_time=1;
while (in*dt)<avg_time
    imu_data=fread(F_IMU,8,'double');
    true_val = fread(F_TRU, 10, 'double');
    imu_err=fread(F_IMUERR, 7, 'double');
    imu_avg=imu_avg+imu_data(3:8);
    in=in+1;
end
imu_avg=imu_avg/in;
%Compute the roll/pitch
%[Cnb E]=align_wander(imu_avg(1:3), [0;0;-1]);
[Cnb E]=alingc_pl_v000(imu_avg(1:3), 0);
Cbn=Cnb';

%Use constant scale for all observations
[Rn, Re, g, sL, cL, WIE_E]=geoparam_v000(pos_g);
Sn=Rn+pos_g(3);
Se=(Re+pos_g(3))*cos(pos_g(1));
S=[Sn 0 0;0 Se 0;0 0 1];

%alignment mode (0=coarse, 1=fine)
%%"mode 2" is standard small angle (PHI-model) model. (I used it to verify the results.)
mode=0;

%Initial covariance
if (mode==0)
    nst=17;
    nst_nav=11;
    P=zeros(nst);
    P(10:11,10:11)=IniErrDef.wander_sP*IniErrDef.wander_sP';
    wander=[sin(0*pi/180);cos(0*pi/180)];
elseif (mode==1)
    nst=16;
    nst_nav=10;
    P=zeros(nst);
    P(10,10)=IniErrDef.att_sP(3,3)*IniErrDef.att_sP(3,3);
    %wander=[sin(105*pi/180);cos(105*pi/180)]; (use the wander provided by
    %the syspath)
    %P(10:11,10:11)=IniErrDef.wander_sP(1,1)+IniErrDef.wander_sP(2,2);
end
P(1:3,1:3)=S*IniErrDef.pos_sP*IniErrDef.pos_sP'*S';
P(4:6,4:6)=IniErrDef.vel_sP*IniErrDef.vel_sP';
P(7:8,7:8)=E(1:2,1:3)*Rimu(1:3,1:3)/avg_time*E(1:2,1:3)';
P(nst_nav+1:nst,nst_nav+1:nst)=sP0imu*sP0imu';
%Cross covarince between acc error and initial attitude
P(7:8,nst_nav+1:nst_nav+3)=E(1:2,1:3)*P(nst_nav+1:nst_nav+3,nst_nav+1:nst_nav+3);
P(nst_nav+1:nst_nav+3,7:8)=P(7:8,nst_nav+1:nst_nav+3)';

%accumulators and control flags
prp_coun=0;
STM_upd=eye(nst);
I=eye(nst);
imu_avg=zeros(6,1);

%%read the input data from the files
imu_data=fread(F_IMU,8,'double');
true_val = fread(F_TRU, 10, 'double');
imu_err=fread(F_IMUERR, 7, 'double');
obs_pos=fread(F_OBS_POS, 4, 'double');

%%Start iterations
while (~feof(F_IMU))
        imu=imu_data(3:8)-ximu;
        pr_coun=imu_data(1);
        [Cbn, vel_n, pos_g, wander]=strapdown_wander_dcm_v000(Cbn, vel_n, pos_g, wander, imu(1:3), imu(4:6), dt);
        imu_avg=imu_avg+imu;
        prp_coun=prp_coun+1;

        if ((~isempty(obs_pos) && pr_coun==obs_pos(1))|| (mod(prp_coun,prp_per)==0)) %%(obs and/or covariance propagation)
            %%First propagate the covariance
            imu_avg=imu_avg/prp_coun;

            %%Form the discrete time navigation model
            if (mode==0 || mode==1)
                [STM Q]=sys_wander_largeheading_v000(pos_g, vel_n, Cbn, wander, imu_avg(1:3), prp_coun*dt, mode, Aimu, Qimu, Cimu, Rimu);
            elseif (mode==2)
                [STM Q]=sys_wander_dcm_v000(pos_g, vel_n, Cbn, wander, imu_avg(1:3), prp_coun*dt, Aimu, Qimu, Cimu, Rimu);
            end

            %Propagate
            P=STM*P*STM'+Q;
            prp_coun=0;
            imu_avg=zeros(6,1);

            %update IMU error states
            ximu=STM(nst_nav+1:nst,nst_nav+1:nst)*ximu;

            if (~isempty(obs_pos)) && (obs_pos(1)==pr_coun) %We have an observation to be processed

                if (mode==0 || mode==1) %switch from coarse/fine to small angle. The new wander angle will be fixed to the last wander estimate.
                    if (mode==0)
                        vda=P(10,10)+P(11,11);
                    elseif (mode==1)
                        vda=P(10,10);
                    end
                    if (vda<(0*pi/180)^2)
                        Cgn=[wander(2) wander(1) 0;-wander(1) wander(2) 0;0 0 1];
                        vel_g=Cgn'*vel_n;
                        mx_a=eye(17);
                        mx_a(10:11,:)=[];
                        mx_a(4:6,10:11)=-[vel_g(2) vel_g(1);-vel_g(1) vel_g(2);0 0];
                        mx_a(9,10:11)=[-wander(2) -wander(1)];

                        if (mode==1)
                            mx_a(:,10)=wander(2)*mx_a(:,10)-wander(1)*mx_a(:,11);
                            mx_a(:,11)=[];
                        end

                        P=mx_a*P*mx_a';

                        mode=2;
                        nst_nav=9;
                        nst=15;
                        I=eye(nst);
                    end
                end

                H=zeros(3,nst); %position update
                H(:,1:3)=eye(3);
                R=S*ObsErrDef.sR*ObsErrDef.sR'*S';
                inno=S*(pos_g-obs_pos(2:4));  %position obs in m

                Rinno_inv=inv(H*P*H'+R);
                K=P*H'*Rinno_inv;

                %correct current states
                dx=K*inno;
                %K(10:11,:)=0;
                %K(12:17,:)=0;
                pos_g=pos_g-S\dx(1:3);
                vel_n=vel_n-dx(4:6);
                mx_a=euler2dcm_v000(dx(7:9));
                Cbn=mx_a*Cbn;
                if (mode==0)    %coarse alignment
                    wander=wander-dx(10:11);
                    %%normalize wander
                    wander=wander/norm(wander);
                elseif (mode==1) %fine alignment
                    wander=[-wander(2) wander(1);wander(1) wander(2)]*[sin(dx(10));cos(dx(10))];
                end
                ximu=ximu+dx(nst_nav+1:nst);
                %update cov
                P=(I-K*H)*P*(I-K*H)'+K*R*K';

                %prepare for the next cycle
                obs_pos=fread(F_OBS_POS, 4, 'double');

                %Check to switch mode
                if (mode==0)    %switch from coarse to fine
                    vda=P(10,10)+P(11,11);
                    if (vda<(0*pi/180)^2)
                        mode=1;
                        nst_nav=10;
                        nst=16;
                        I=eye(nst);

                        mx_a=eye(17,17);
                        mx_a(11,:)=[];
                        mx_a(10,10)=wander(2);
                        mx_a(10,11)=-wander(1);
                        P=mx_a*P*mx_a';
                    end
                end

            end

            %%Write the results (Cbg, Vb, pos_g)
            vr_a=diag(P).^0.5;

            imu_err=imu_err(2:7)-ximu;
            imu_std=vr_a(nst_nav+1:end);

            pos_err=S*(pos_g-true_val(2:4));
            pos_std=vr_a(1:3);

            vel_err=Cbn'*vel_n-true_val(5:7);
            mx_a=[Cbn' -Cbn'*skew(vel_n)];
            mx_b=P(4:9,4:9);
            mx_c=mx_a*mx_b*mx_a';
            vel_std=diag(mx_c).^0.5;
%             Cgn=[wander(2) wander(1) 0;-wander(1) wander(2) 0;0 0 1];
%             Cbg_true=euler2dcm_v000(true_val(8:10));
%             vel_err=Cgn'*vel_n-Cbg_true*true_val(5:7);

            Cbg_true=euler2dcm_v000(true_val(8:10));
            Cgn=[wander(2) wander(1) 0;-wander(1) wander(2) 0;0 0 1];
            Cbg=Cgn'*Cbn;
            mx_a=eye(3)-Cbg*Cbg_true';
            att_err=[-mx_a(2,3);mx_a(1,3);-mx_a(1,2)]; %-dcm2euler_v000(Cgn'*Cbn*Cbn_true');
            if (mode==0)
                mx_a=[Cgn' [0 0;0 0;-wander(2) wander(1)]];
            elseif (mode==1)
                mx_a=[Cgn' -[0;0;1]];
            elseif (mode==2)
                mx_a=Cgn';
            end
            mx_b=P(7:nst_nav,7:nst_nav);
            mx_c=mx_a*mx_b*mx_a';
            att_std=diag(mx_c).^0.5;

%             %%%Write the results (Cbn, Vn, pos_g) (only for mode 0 and
%             %%%mode1 - revise for mode2)
%             vr_a=diag(P).^0.5;
%
%             imu_err=imu_err(2:7)-ximu;
%             imu_std=vr_a(nst_nav+1:end);
%
%             pos_err=S*(pos_g-true_val(2:4));
%             pos_std=vr_a(1:3);
%
%             Cbg_true=euler2dcm_v000(true_val(8:10));
%             Cgn=[wander(2) wander(1) 0;-wander(1) wander(2) 0;0 0 1];
%             Cbn_true=Cgn*Cbg_true;
%             vel_err=vel_n-Cbn_true*true_val(5:7);
%             if (mode==0)
%                 vr_a=Cbg_true*true_val(5:7);
%                 mx_a=[eye(3);[vr_a(2) vr_a(1);-vr_a(1) vr_a(2);0 0]];
%                 mx_b=P([4 5 6 10 11],[4 5 6 10 11]);
%             elseif (mode==1)
%                 vr_a=Cbn_true*true_val(5:7);
%                 mx_a=[eye(3);[vr_a(2);-vr_a(1);0]];
%                 mx_b=P([4 5 6 10],[4 5 6 10]);
%             end
%             mx_c=mx_a*mx_b*mx_a';
%             vel_std=diag(mx_c).^0.5;
%
%
%             mx_a=eye(3)-Cbn*Cbn_true';
%             att_err=[-mx_a(2,3);mx_a(1,3);-mx_a(1,2)]; %-dcm2euler_v000(Cgn'*Cbn*Cbn_true');
%             if (mode==0)
%                 mx_a=[eye(3) [0 0;0 0;-wander(2) -wander(1)]];
%             elseif (mode==1)
%                 mx_a=[eye(3) -[0;0;1]];
%             end
%             mx_b=P(7:nst_nav,7:nst_nav);
%             mx_c=mx_a*mx_b*mx_a';
%             att_std=diag(mx_c).^0.5;

            fwrite(F_ERR,[pr_coun;pos_err;vel_err;att_err;imu_err;pos_std;vel_std;att_std;imu_std],'double');
            fwrite(F_RES, [pr_coun;pos_g; Cbn'*vel_n; dcm2euler_v000(Cgn'*Cbn);atan2(wander(1),wander(2));ximu],'double');
        end

        %Read the next values
        true_val = fread(F_TRU, 10, 'double');
        imu_err=fread(F_IMUERR, 7, 'double');
        imu_data=fread(F_IMU,8,'double');
end

fclose all;
err=readbin_v000([PATH 'err.bin'],31);
%imu_err=readbin_v000([PATH '\imuerr.bin'],7);
res=readbin_v000([PATH 'res.bin'],17);
tru=readbin_v000([PATH 'mnav.bin'],10);

% return;
% figure;
% plot(res(2,:), res(3,:))
% hold on
% plot(tru(2,:),(tru(3,:)),'r--');
% grid;

figure;
i=9;
plot(abs(err(1+i,:))*180/pi);
hold on;
plot(err(16+i,:)*180/pi,'r');
grid;

% figure;
% i=4;
% plot(abs(err(1+i,:)));
% hold on;
% plot(err(16+i,:),'r');
% grid;

repo: ekfmonoslam


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

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

% Last Modified by GUIDE v2.5 13-Feb-2013 04:22:47

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

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

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

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes FUN_PLOTTER wait for user response (see UIRESUME)
% uiwait(handles.figure1);

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

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global z thta
cla(handles.axes1)
reset(handles.axes1)
set(handles.axes1,'box','on')
set(handles.axes1,'xtick',[])
set(handles.axes1,'ytick',[])
z = 'sin(x)';
v = get(handles.checkbox2,'value');
if v == 1
    thta = 0:3*360;
    axes(handles.axes1)
    plot(thta,sind(thta),'linewidth',3)
    set(handles.text1,'string','SINE PLOT FROM 0 TO 6*PI')
    grid on
else
    msgbox('PLZ SET RANGE...TRY AGAIN','WARNING','warn','modal')
    set(handles.text1,'string','')
    cla(handles.axes1)
    reset(handles.axes1)
    set(handles.axes1,'box','on')
    set(handles.axes1,'xtick',[])
    set(handles.axes1,'ytick',[])
    return
end

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global z thta
cla(handles.axes1)
reset(handles.axes1)
set(handles.axes1,'box','on')
set(handles.axes1,'xtick',[])
set(handles.axes1,'ytick',[])
z = 'cos(x)';
v = get(handles.checkbox2,'value');
if v == 1
    thta = 0:3*360;
    axes(handles.axes1)
    plot(thta,cosd(thta),'linewidth',3)
    set(handles.text1,'string','COSINE PLOT FROM 0 TO 6*PI')
    grid on
else
    msgbox('PLZ SET RANGE...TRY AGAIN','WARNING','warn','modal')
    set(handles.text1,'string','')
    cla(handles.axes1)
    reset(handles.axes1)
    set(handles.axes1,'box','on')
    set(handles.axes1,'xtick',[])
    set(handles.axes1,'ytick',[])
    return
end

% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global z thta
cla(handles.axes1)
reset(handles.axes1)
set(handles.axes1,'box','on')
set(handles.axes1,'xtick',[])
set(handles.axes1,'ytick',[])
z = 'tan(x)';
v = get(handles.checkbox2,'value');
if v == 1
    thta = 0:3*360;
    axes(handles.axes1)
    plot(thta,tand(thta),'linewidth',3)
    set(handles.text1,'string','TANGENT PLOT FROM 0 TO 6*PI')
    grid on
else
    msgbox('PLZ SET RANGE...TRY AGAIN','WARNING','warn','modal')
    set(handles.text1,'string','')
    cla(handles.axes1)
    reset(handles.axes1)
    set(handles.axes1,'box','on')
    set(handles.axes1,'xtick',[])
    set(handles.axes1,'ytick',[])
    return
end

% --- Executes on button press in pushbutton4.
function pushbutton4_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global z thta
cla(handles.axes1)
reset(handles.axes1)
set(handles.axes1,'box','on')
set(handles.axes1,'xtick',[])
set(handles.axes1,'ytick',[])
z = 'cot(x)';
v = get(handles.checkbox2,'value');
if v == 1
    thta = 0:3*360;
    axes(handles.axes1)
    plot(thta,cotd(thta),'linewidth',3)
    set(handles.text1,'string','COTANGENT PLOT FROM 0 TO 6*PI')
    grid on
else
    msgbox('PLZ SET RANGE...TRY AGAIN','WARNING','warn','modal')
    set(handles.text1,'string','')
    cla(handles.axes1)
    reset(handles.axes1)
    set(handles.axes1,'box','on')
    set(handles.axes1,'xtick',[])
    set(handles.axes1,'ytick',[])
    return
end

% --- Executes on button press in pushbutton5.
function pushbutton5_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global z thta
cla(handles.axes1)
reset(handles.axes1)
set(handles.axes1,'box','on')
set(handles.axes1,'xtick',[])
set(handles.axes1,'ytick',[])
z = 'sec(x)';
v = get(handles.checkbox2,'value');
if v == 1
    thta = 0:3*360;
    axes(handles.axes1)
    plot(thta,secd(thta),'linewidth',3)
    set(handles.text1,'string','SECANT PLOT FROM 0 TO 6*PI')
    grid on
else
    msgbox('PLZ SET RANGE...TRY AGAIN','WARNING','warn','modal')
    set(handles.text1,'string','')
    cla(handles.axes1)
    reset(handles.axes1)
    set(handles.axes1,'box','on')
    set(handles.axes1,'xtick',[])
    set(handles.axes1,'ytick',[])
    return
end

% --- Executes on button press in pushbutton6.
function pushbutton6_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global z thta
cla(handles.axes1)
reset(handles.axes1)
set(handles.axes1,'box','on')
set(handles.axes1,'xtick',[])
set(handles.axes1,'ytick',[])
z = 'csc(x)';
v = get(handles.checkbox2,'value');
if v == 1
    thta = 0:3*360;
    axes(handles.axes1)
    plot(thta,cscd(thta),'linewidth',3)
    set(handles.text1,'string','COSECANT PLOT FROM 0 TO 6*PI')
    grid on
else
    msgbox('PLZ SET RANGE...TRY AGAIN','WARNING','warn','modal')
    set(handles.text1,'string','')
    cla(handles.axes1)
    reset(handles.axes1)
    set(handles.axes1,'box','on')
    set(handles.axes1,'xtick',[])
    set(handles.axes1,'ytick',[])
    return
end

% --- Executes on button press in pushbutton7.
function pushbutton7_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton7 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global z x
cla(handles.axes1)
reset(handles.axes1)
set(handles.axes1,'box','on')
set(handles.axes1,'xtick',[])
set(handles.axes1,'ytick',[])
z = 'log(t)';
v = get(handles.checkbox2,'value');
if v == 1
    x = 0:0.25:10;
    axes(handles.axes1)
    plot(x,log(x),'linewidth',3)
    set(handles.text1,'string','LOGRITHMIC PLOT FROM 0 TO 10')
    grid on
else
    msgbox('PLZ SET RANGE...TRY AGAIN','WARNING','warn','modal')
    set(handles.text1,'string','')
    cla(handles.axes1)
    reset(handles.axes1)
    set(handles.axes1,'box','on')
    set(handles.axes1,'xtick',[])
    set(handles.axes1,'ytick',[])
    return
end

% --- Executes on button press in pushbutton8.
function pushbutton8_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton8 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global z x
cla(handles.axes1)
reset(handles.axes1)
set(handles.axes1,'box','on')
set(handles.axes1,'xtick',[])
set(handles.axes1,'ytick',[])
z = 'exp(t)';
v = get(handles.checkbox2,'value');
if v == 1
    x = 0:0.25:10;
    axes(handles.axes1)
    plot(x,exp(x),'linewidth',3)
    set(handles.text1,'string','EXPONENTIAL PLOT FROM 0 TO 10')
    grid on
else
    msgbox('PLZ SET RANGE...TRY AGAIN','WARNING','warn','modal')
    set(handles.text1,'string','')
    cla(handles.axes1)
    reset(handles.axes1)
    set(handles.axes1,'box','on')
    set(handles.axes1,'xtick',[])
    set(handles.axes1,'ytick',[])
    return
end

% --- Executes on button press in checkbox2.
function checkbox2_Callback(hObject, eventdata, handles)
% hObject    handle to checkbox2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of checkbox2
set(handles.checkbox3,'value',0)
set(handles.edit1,'enable','inactive')
set(handles.edit2,'enable','inactive')
% --- Executes on button press in checkbox3.
function checkbox3_Callback(hObject, eventdata, handles)
% hObject    handle to checkbox3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of checkbox3
set(handles.checkbox2,'value',0)
set(handles.edit1,'enable','on')
set(handles.edit2,'enable','on')

function edit1_Callback(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit1 as text
%        str2double(get(hObject,'String')) returns contents of edit1 as a double

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

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

function edit2_Callback(hObject, eventdata, handles)
% hObject    handle to edit2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit2 as text
%        str2double(get(hObject,'String')) returns contents of edit2 as a double

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

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

function edit3_Callback(hObject, eventdata, handles)
% hObject    handle to edit3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit3 as text
%        str2double(get(hObject,'String')) returns contents of edit3 as a double

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

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

% --- Executes on button press in pushbutton10.
function pushbutton10_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton10 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global z p1
v = get(handles.checkbox2,'value');
if v == 1
    msgbox('CAN NOT WORK ON DEFAULT RANGE','WARNING','warn','modal')
    return
end
v1 = get(handles.checkbox3,'value');
if v1 == 0
    msgbox('FIRST SET RANGE','WARNING','warn','modal')
    return
end
a = str2double(get(handles.edit1,'string'));
b = str2double(get(handles.edit2,'string'));
if a > b
    msgbox('MIN IS GR8ER THAN MAX','WARNING','warn','modal')
    set(handles.edit1,'string','')
    set(handles.edit2,'string','')
    return
end
fx2 = get(handles.edit3,'string');
z = fx2;
if isempty(fx2) || isempty(a) || isempty(b)
    msgbox('EITHER MAX OR MIN OR FUNCTION FIELD IS EMPTY OR INVALID','WARNING...','modal')
    set(handles.edit1,'string','')
    set(handles.edit2,'string','')
    set(handles.edit3,'string','')
    return
end
fx3 = inline(fx2);
p1 = a:0.25:b;
p2 = p1;
for l = 1:length(p2)
    mag(l) = fx3(p2(l));
end
axes(handles.axes1)
plot(p2,mag,'linewidth',3)
grid on
set(handles.text1,'string',strcat('PLOT OF ',fx2,' FROM RANGE ',num2str(a),' TO ',num2str(b)))

% --- Executes on button press in pushbutton11.
function pushbutton11_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton11 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global z thta x p1
d = z;
k = z;
switch k
    case 'sin(x)'
        fx = sym(k);
        n = str2double(inputdlg('INPUT ORDER OF DIFFRENTIATION','DY/DX'));
        if isnan(n)
            return
        end
        fx1 = diff(fx,n);
        for l = 1:length(thta)
            mag(l) = subs(fx1,deg2rad(thta(l)));
        end
        axes(handles.axes1)
        hold on
        plot(thta,mag,'r','linewidth',3)
        grid on
        %legend('sin\theta','d(sin\theta)/dx')
        set(handles.text1,'string','DIFFRENTIATION OF SIN(X) AND SIN(X)')
    case 'cos(x)'
        fx = sym(k);
        n = str2double(inputdlg('INPUT ORDER OF DIFFRENTIATION','DY/DX'));
        if isnan(n)
            return
        end
        fx1 = diff(fx,n);
        for l = 1:length(thta)
            mag(l) = subs(fx1,deg2rad(thta(l)));
        end
        axes(handles.axes1)
        hold on
        plot(thta,mag,'r','linewidth',3)
        grid on
        %legend('cos\theta','d(cos\theta)/dx')
        set(handles.text1,'string','DIFFRENTIATION OF COS(X) AND COS(X)')
    case 'tan(x)'
        fx = sym(k);
        fx1 = diff(fx);
        for l = 1:length(thta)
            mag(l) = subs(fx1,deg2rad(thta(l)));
        end
        axes(handles.axes1)
        hold on
        plot(thta,mag,'r','linewidth',3)
        grid on
        %legend('tan\theta','d(tan\theta)/dx')
        set(handles.text1,'string','DIFFRENTIATION OF TAN(X) AND TAN(X)')
    case 'cot(x)'
        fx = sym(k);
        fx1 = diff(fx);
        for l = 1:length(thta)
            mag(l) = subs(fx1,deg2rad(thta(l)));
        end
        axes(handles.axes1)
        hold on
        plot(thta,mag,'r','linewidth',3)
        grid on
        %legend('cot\theta','d(cot\theta)/dx')
        set(handles.text1,'string','DIFFRENTIATION OF COT(X) AND COT(X)')
    case 'sec(x)'
        fx = sym(k);
        fx1 = diff(fx);
        for l = 1:length(thta)
            mag(l) = subs(fx1,deg2rad(thta(l)));
        end
        axes(handles.axes1)
        hold on
        plot(thta,mag,'r','linewidth',3)
        grid on
        %legend('sec\theta','d(sec\theta)/dx')
        set(handles.text1,'string','DIFFRENTIATION OF SEC(X) AND SEC(X)')
    case 'csc(x)'
        fx = sym(k);
        fx1 = diff(fx);
        for l = 1:length(thta)
            mag(l) = subs(fx1,deg2rad(thta(l)));
        end
        axes(handles.axes1)
        hold on
        plot(thta,mag,'r','linewidth',3)
        grid on
        %legend('csc\theta','d(csc\theta)/dx')
        set(handles.text1,'string','DIFFRENTIATION OF CSC(X) AND CSC(X)')
    case 'log(t)'
        fx = sym(k);
        fx1 = diff(fx);
        r = x;
        for l = 1:length(r)
            mag(l) = subs(fx1,r(l));
        end
        axes(handles.axes1)
        hold on
        plot(r,mag,'r','linewidth',3)
        grid on
        %legend('log(x)','d(log(x))/dx')
        set(handles.text1,'string','DIFFRENTIATION OF LOG(X) AND LOG(X)')
    case 'exp(t)'
        fx = sym(k);
        fx1 = diff(fx);
        r = x;
        for l = 1:length(r)
            mag(l) = subs(fx1,r(l));
        end
        axes(handles.axes1)
        hold on
        plot(r,mag,'r','linewidth',1.5)
        grid on
        %legend('e^x','d(e^x)/dx')
        set(handles.text1,'string','DIFFRENTIATION OF EXP(X) AND EXP(X)')
    case d
        fx = sym(k);
        fx1 = diff(fx);
        r = p1;
        for l = 1:length(r)
            mag(l) = subs(fx1,r(l));
        end
        axes(handles.axes1)
        hold on
        plot(r,mag,'r','linewidth',1.8)
        grid on
        %legend(d,strcat('d(',d,')/dx'))
        set(handles.text1,'string',strcat('DIFFRENTIATION OF ',d,' AND ',d))
end

% --- Executes on button press in pushbutton12.
function pushbutton12_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton12 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global z thta x p1
d = z;
k = z;
switch k
    case 'sin(x)'
        fx = sym(k);
        fx1 = int(fx);
        for l = 1:length(thta)
            mag(l) = subs(fx1,deg2rad(thta(l)));
        end
        axes(handles.axes1)
        hold on
        plot(thta,mag,'g','linewidth',3)
        grid on
        %legend('sin\theta','d(sin\theta)/dx')
        set(handles.text1,'string','INTEGRATION OF SIN(X) AND SIN(X)')
    case 'cos(x)'
        fx = sym(k);
        fx1 = int(fx);
        for l = 1:length(thta)
            mag(l) = subs(fx1,deg2rad(thta(l)));
        end
        axes(handles.axes1)
        hold on
        plot(thta,mag,'g','linewidth',3)
        grid on
        %legend('cos\theta','d(cos\theta)/dx')
        set(handles.text1,'string','INTEGRATION OF COS(X) AND COS(X)')
    case 'tan(x)'
        fx = sym(k);
        fx1 = int(fx);
        for l = 1:length(thta)
            mag(l) = subs(fx1,deg2rad(thta(l)));
        end
        axes(handles.axes1)
        hold on
        plot(thta,mag,'g','linewidth',3)
        grid on
        %legend('tan\theta','d(tan\theta)/dx')
        set(handles.text1,'string','INTEGRATION OF TAN(X) AND TAN(X)')
    case 'cot(x)'
        fx = sym(k);
        fx1 = int(fx);
        for l = 1:length(thta)
            mag(l) = subs(fx1,deg2rad(thta(l)));
        end
        axes(handles.axes1)
        hold on
        plot(thta,mag,'g','linewidth',3)
        grid on
        %legend('cot\theta','d(cot\theta)/dx')
        set(handles.text1,'string','INTEGRATION OF COT(X) AND COT(X)')
    case 'sec(x)'
        fx = sym(k);
        fx1 = int(fx);
        for l = 1:length(thta)
            mag(l) = subs(fx1,deg2rad(thta(l)));
        end
        axes(handles.axes1)
        hold on
        plot(thta,mag,'g','linewidth',3)
        grid on
        %legend('sec\theta','d(sec\theta)/dx')
        set(handles.text1,'string','INTEGRATION OF SEC(X) AND SEC(X)')
    case 'csc(x)'
        fx = sym(k);
        fx1 = int(fx);
        for l = 1:length(thta)
            mag(l) = subs(fx1,deg2rad(thta(l)));
        end
        axes(handles.axes1)
        hold on
        plot(thta,mag,'g','linewidth',3)
        grid on
        %legend('csc\theta','d(csc\theta)/dx')
        set(handles.text1,'string','INTEGRATION OF CSC(X) AND CSC(X)')
    case 'log(t)'
        fx = sym(k);
        fx1 = int(fx);
        r = x;
        for l = 1:length(r)
            mag(l) = subs(fx1,r(l));
        end
        axes(handles.axes1)
        hold on
        plot(r,mag,'g','linewidth',3)
        grid on
        %legend('log(x)','d(log(x))/dx')
        set(handles.text1,'string','INTEGRATION OF LOG(X) AND LOG(X)')
    case 'exp(t)'
        fx = sym(k);
        fx1 = int(fx);
        r = x;
        for l = 1:length(r)
            mag(l) = subs(fx1,r(l));
        end
        axes(handles.axes1)
        hold on
        plot(r,mag,'g','linewidth',1.5)
        grid on
        %legend('e^x','d(e^x)/dx')
        set(handles.text1,'string','INTEGRATION OF EXP(X) AND EXP(X)')

     case d
        fx = sym(k);
        fx1 = int(fx);
        r = p1;
        for l = 1:length(r)
            mag(l) = subs(fx1,r(l));
        end
        axes(handles.axes1)
        hold on
        plot(r,mag,'g','linewidth',2.3)
        grid on
        %legend(d,strcat('d(',d,')/dx'))
        set(handles.text1,'string',strcat('INTEGRATION OF ',d,' AND ',d))
end

% --------------------------------------------------------------------
function CL_EAR_Callback(hObject, eventdata, handles)
% hObject    handle to CL_EAR (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
cla(handles.axes1)
reset(handles.axes1)
set(handles.axes1,'box','on')
set(handles.axes1,'xtick',[])
set(handles.axes1,'ytick',[])
set(handles.edit1,'string','')
set(handles.edit2,'string','')
set(handles.edit3,'string','')
set(handles.checkbox2,'value',0)
set(handles.checkbox3,'value',0)

% --------------------------------------------------------------------
function EXI_T_Callback(hObject, eventdata, handles)
% hObject    handle to EXI_T (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
close

% --------------------------------------------------------------------
function CLEAR_AXES_Callback(hObject, eventdata, handles)
% hObject    handle to CLEAR_AXES (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
cla(handles.axes1)
reset(handles.axes1)
set(handles.axes1,'box','on')
set(handles.axes1,'xtick',[])
set(handles.axes1,'ytick',[])


text:
%% PRIME_SCRIPT is a script to call PRIME_FUN.
%
%  Discussion:
%
%    The BATCH command runs scripts, not functions.  So we have to write
%    this short script if we want to work with BATCH!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2012
%
%  Author:
%
%    John Burkardt
%
  n = 10000000;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PRIME_SCRIPT\n' );
  fprintf ( 1, '  Call PRIME_FUN to count prime numbers from 1 to %d\n', n );

  tic
  total = prime_fun ( n );
  toc


text:
function v = getOr(s, field, default)
%getOr Returns the structure field or a default if either don't exist
%   v = getOr(s, field, [default]) returns the 'field' of the structure 's'
%   or 'default' if the structure is empty of the field does not exist. If
%   default is not specified it defaults to []. 'field' can also be a cell
%   array of fields, in which case it will check for all of them and return
%   the value of the first field that exists, if any (otherwise the default
%   value).

if nargin < 3
  default = [];
end

fieldExists = isfield(s, field);
if any(fieldExists)
  if iscellstr(field)
    v = s.(field{find(fieldExists, 1)});
  else
    v = s.(field);
  end
else
  v = default;
end

end

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
classdef DPCA < Algorithms.CDAlg
    properties
        nPCUsed = uint8(255);
    end
    methods
        function obj = DPCA(k)
            obj@Algorithms.CDAlg('Differential Principal Component Analysis', '');
            if exist('k', 'var'), obj.nPCUsed=uint8(k); end
        end
        [DI, k] = detectChange(obj, t1, t2);
    end
end

author: Bobholamovic

repo: ChangeDetectionToolbox


text:
% realproba() - compute the effective probability of the value
%               in the sample.
%
% Usage:
%   >> [probaMap, probaDist ] = realproba( data, discret);
%
% Inputs:
%   data       - the data onto which compute the probability
%   discret    - discretisation factor (default: (size of data)/5)
%                if 0 base the computation on a Gaussian
%                approximation of the data
%
% Outputs:
%   probaMap   - the probabilities associated with the values
%   probaDist  - the probabilities distribution
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001

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

function [ probaMap, sortbox ] = realproba( data, bins );

if nargin < 1
	help realproba;
	return;
end;
if nargin < 2
	bins = round(size(data,1)*size(data,2)/5);
end;

if bins > 0
	% COMPUTE THE DENSITY FUNCTION
	% ----------------------------
	SIZE = size(data,1)*size(data,2);
	sortbox = zeros(1,bins);
	minimum =  min(data(:));
	maximum =  max(data(:));
	data = floor((data - minimum )/(maximum - minimum)*(bins-1))+1;
    if any(any(isnan(data))), warning('Binning failed - could be due to zeroed out channel'); end;
	for index=1:SIZE
		sortbox(data(index)) = sortbox(data(index))+1;
	end;
	probaMap = sortbox(data) / SIZE;
	sortbox  = sortbox / SIZE;
else
	% BASE OVER ERROR FUNCTION
	% ------------------------
	data     = (data-mean(data(:)))./std(data(:));
	probaMap = exp(-0.5*( data.*data ))/(2*pi);
	probaMap = probaMap/sum(probaMap); % because the total surface under a normalized Gaussian is 2
	sortbox  = probaMap/sum(probaMap);
end;
return;

author: PatternRecognition

repo: OpenBMI


text:
function S1 = my_conv2(S1, sig, varargin)
% S1 is the matrix to be filtered along a choice of axes
% sig is either a scalar or a sequence of scalars, one for each axis to be filtered
% varargin can be the dimensions to do filtering, if len(sig) != x.shape
% if sig is scalar and no axes are provided, the default axis is 2

if sig>.25
    idims = 2;
    if ~isempty(varargin)
        idims = varargin{1};
    end
    if numel(idims)>1 && numel(sig)>1
        sigall = sig;
    else
        sigall = repmat(sig, numel(idims), 1);
    end

    for i = 1:length(idims)
        sig = sigall(i);

        idim = idims(i);
        Nd = ndims(S1);

        S1 = permute(S1, [idim 1:idim-1 idim+1:Nd]);

        dsnew = size(S1);

        S1 = reshape(S1, size(S1,1), []);
        dsnew2 = size(S1);

        % NN = size(S1,1);
        % NT = size(S1,2);

        tmax = ceil(4*sig);
        dt = -tmax:1:tmax;
        gaus = exp( - dt.^2/(2*sig^2));
        gaus = gaus'/sum(gaus);

        % Norms = conv(ones(NT,1), gaus, 'same');
        % Smooth = zeros(NN, NT);
        % for n = 1:NN
        %    Smooth(n,:) = (conv(S1(n,:)', gaus, 'same')./Norms)';
        % end

        cNorm = filter(gaus, 1, cat(1, ones(dsnew2(1), 1), zeros(tmax,1)));
        cNorm = cNorm(1+tmax:end, :);
        S1 = filter(gaus, 1, cat(1, S1, zeros([tmax, dsnew2(2)])));
        S1(1:tmax, :) = [];
        S1 = reshape(S1, dsnew);

        S1 = bsxfun(@rdivide, S1, cNorm);

        S1 = permute(S1, [2:idim 1 idim+1:Nd]);
    end
end