text:
function localangle = point_angle(bw, seed, R, start, NeighborNum)

% This function computes the angle of seed from the start point
if (nargin == 4)
    NeighborNum =3;
end

[M, N]= size(bw);

% image(y,x) <==>(x-1)*M + y
seedy = mod(seed, M);
if (seedy==0) seedy = M; end
seedx = 1 + (seed - seedy)/M;
starty = mod(start, M);
if (starty==0) starty = M; end
startx = 1 + (start - starty)/M;

dy = 1 + R + starty - seedy;
dx = 1 + R + startx - seedx;

% set the region border from 1:8*R
region = zeros(2*R+1);
region(R+1:-1:1,end) = (1:R+1)';
region(1,end-1:-1:2)= R+2:3*R;
region(1:end,1)= (3*R+1:5*R+1)';
region(end,2:end-1) = 5*R+2:7*R;
region(end:-1:R+2,end) = (7*R+1:8*R)';

startidx = region(dy, dx);
anglevec = point_anglevec(bw, seed, R);

if (startidx==0)
    localangle = [];
    return;
end

if (anglevec(startidx)==0)
    localangle = [];
    return;
end

localangle = findangle(anglevec, startidx);
localnum = prod(size(localangle));

if (localnum > NeighborNum)
    localangle = localangle(1:NeighborNum);
end
if (localnum < NeighborNum)
    localangle = [localangle, zeros(1, NeighborNum - localnum )];
end

author: Sable


text:
function disp( x, prefix ) %#ok
if nargin < 2, prefix = ''; end
disp( [ prefix, 'cvx base object' ] );

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function x = cs_usolve (U,b)                                                %#ok
%CS_USOLVE solve a sparse upper triangular system U*x=b.
%   x = cs_usolve(U,b) computes x = U\b, U must be lower triangular with a
%   zero-free diagonal.  b must be a column vector.  x is full if b is full.
%   If b is sparse, x is sparse but nonzero pattern of x is NOT sorted (it is
%   returned in topological order).
%
%   Example:
%       Prob = UFget ('HB/arc130') ; A = Prob.A ; n = size (A,1) ;
%       b = rand (n,1);
%       [L U p q] = cs_lu (A) ;
%       x = cs_usolve (U, cs_lsolve (L, b(p))) ;   % x = U \ (L \ b(p)) ;
%       x (q) = x ;
%       norm (A*x-b)
%
%   See also CS_LSOLVE, CS_LTSOLVE, CS_UTSOLVE, MLDIVIDE.

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

error ('cs_usolve mexFunction not found') ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function dcline = psse_convert_hvdc(dc, bus)
%PSSE_CONVERT_HVDC Convert HVDC data from PSS/E RAW to MATPOWER
%   DCLINE = PSSE_CONVERT_HVDC(DC, BUS)
%
%   Convert all two terminal HVDC line data read from a PSS/E
%   RAW data file into MATPOWER format. Returns a dcline matrix for
%   inclusion in a MATPOWER case struct.
%
%   Inputs:
%       DC  : matrix of raw two terminal HVDC line data returned by
%             PSSE_READ in data.twodc.num
%       BUS : MATPOWER bus matrix
%
%   Output:
%       DCLINE : a MATPOWER dcline matrix suitable for inclusion in
%                a MATPOWER case struct.
%
%   See also PSSE_CONVERT.

%   MATPOWER
%   Copyright (c) 2014-2016, Power Systems Engineering Research Center (PSERC)
%   by Yujia Zhu, PSERC ASU
%   and Ray Zimmerman, PSERC Cornell
%   Based on mpdcin.m and mpqhvdccal.m, written by:
%       Yujia Zhu, Jan 2014, yzhu54@asu.edu.
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% define named indices into bus, gen, branch matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
c = idx_dcline;

nb = size(bus, 1);
ndc = size(dc, 1);
e2i = sparse(bus(:, BUS_I), ones(nb, 1), 1:nb, max(bus(:, BUS_I)), 1);
if ~ndc
    dcline = [];
    return;
end

%% extract data
MDC = dc(:,2); % Control mode
SETVL = dc(:,4); % depend on control mode: current or power demand
VSCHD = dc(:,5); % scheduled compounded dc voltage
ANMXR = dc(:,15); % nominal maximum rectifier firing angle
ANMNR = dc(:,16); % nominal minimum rectifier firing angle
GAMMX = dc(:,32); % nominal maximum inverter firing angle
GAMMN = dc(:,33); % nominal minimum inverter firing angle
SETVL = abs(SETVL);
% Convert the voltage on rectifier side and inverter side
% The value is calculated as basekV/VSCHD
% basekV is the bus base voltage, VSCHD is the scheduled compounded
% voltage
dcline = zeros(ndc, c.LOSS1); % initiate the hvdc data format
indr = dc(:,13); % rectifier end bus number
indi = dc(:,30); % inverter end bus number
dcind = [indr indi];
% bus nominal voltage
Vr = bus(e2i(indr), VM);
Vi = bus(e2i(indi), VM);
%% Calculate the real power input at the from end
PMW = zeros(ndc, 1);
for i = 1:ndc
    if MDC(i) == 1
        PMW(i) = SETVL(i); % SETVL is the desired real power demand
    elseif MDC(i) == 2;
        PMW(i) = SETVL(i)*VSCHD(i)/1000; % SETVL is the current in amps (need devide 1000 to convert to MW)
    else PMW(i) = 0;
    end
end
%% calculate reactive power limits
[Qrmin,Qrmax] = psse_convert_hvdc_Qlims(ANMXR,ANMNR,PMW);    %% rectifier end
[Qimin,Qimax] = psse_convert_hvdc_Qlims(GAMMX,GAMMN,PMW);    %% inverter end
%% calculate the loss coefficient (Only consider the l1)
% l1 = P'.*RDC;

%% conclude all info
status = ones(ndc, 1);
status(MDC==0) = 0;     %% set status of blocked HVDC lines to zero
% dcline(:,[1 2 3 4 5 8 9 10 11 12 13 14 15]) = [indr,indi,status,PMW, PMW, Vr, Vi,0.85*PMW, 1.15*PMW, Qrmin, Qrmax, Qimin, Qimax];
dcline(:, [c.F_BUS c.T_BUS c.BR_STATUS c.PF c.PT c.VF c.VT ...
            c.PMIN c.PMAX c.QMINF c.QMAXF c.QMINT c.QMAXT]) = ...
    [indr indi status PMW PMW Vr Vi 0.85*PMW 1.15*PMW Qrmin Qrmax Qimin Qimax];

function [Qmin, Qmax] = psse_convert_hvdc_Qlims(alphamax,alphamin,P)
%PSSE_CONVERT_HVDC_QLIMS calculate HVDC line reactive power limits
%
%   [Qmin, Qmax] = psse_convert_hvdc_Qlims(alphamax,alphamin,P)
%
% Inputs:
%       alphamax :  maximum firing angle
%       alphamin :  minimum steady-state rectifier firing angle
%       P :         real power demand
% Outputs:
%       Qmin :  lower limit of reactive power
%       Qmax :  upper limit of reactive power
%
% Note:
%   This function calculates the reactive power at the rectifier or inverter
%   end. It is assumed the maximum overlap angle is 60 degree (see
%   Kimbark's book). The maximum reactive power is calculated with the
%   power factor:
%       pf = acosd(0.5*(cosd(alphamax(i))+cosd(60))),
%   where, 60 is the maximum delta angle.

len = length(alphamax);
phi = zeros(size(alphamax));
Qmin = phi;
Qmax = phi;
for i = 1:len
    %% minimum reactive power calculated under assumption of no overlap angle
    %% i.e. power factor equals to tan(alpha)
    Qmin(i) = P(i)*tand(alphamin(i));

    %% maximum reactive power calculated when overlap angle reaches max
    %% value (60 deg). I.e.
    %%      cos(phi) = 1/2*(cos(alpha)+cos(delta))
    %%      Q = P*tan(phi)
    phi(i) = acosd(0.5*(cosd(alphamax(i))+cosd(60)));
    Qmax(i) = P(i)*tand(phi(i));
    if Qmin(i)<0
        Qmin(i) = -Qmin(i);
    end
    if Qmax(i)<0
        Qmax(i) = -Qmax(i);
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% op_freqrange.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% out=op_freqrange(in,ppmmin,ppmmax);
%
% DESCRIPTION:
% Output only a specified frequency range of the input spectrum.
%
% INPUTS:
% in         = input data in matlab structure format.
% ppmmin     = minimum extent of frequency range in ppm.
% ppmmax     = maximum extent of frequency range in ppm.
%
% OUTPUTS:
% out        = Output following frequency range selection.

function out=op_freqrange(in,ppmmin,ppmmax);

%Calculate Specs using fft
fullspecs=fftshift(ifft(in.fids,[],in.dims.t),in.dims.t);

%now take only the specified range of the spectrum
specs=fullspecs(in.ppm>ppmmin & in.ppm<ppmmax,:,:);

%convert back to time domain
%if the length of Fids is odd, then you have to do a circshift of one to
%make sure that you don't introduce a small frequency shift into the fids
%vector.
if mod(size(specs,in.dims.t),2)==0
    %disp('Length of vector is even.  Doing normal conversion');
    fids=fft(fftshift(specs,in.dims.t),[],in.dims.t);
else
    %disp('Length of vector is odd.  Doing circshift by 1');
    fids=fft(circshift(fftshift(specs,in.dims.t),1),[],in.dims.t);
end

%calculate the size;
sz=size(fids);

%calculate the ppm scale
ppm=in.ppm(in.ppm>ppmmin & in.ppm<ppmmax);

%calculate the new spectral width and dwelltime:
dppm=abs(ppm(2)-ppm(1));
ppmrange=abs((ppm(end)-ppm(1)))+dppm;
spectralwidth=ppmrange*in.Bo*42.577;
dwelltime=1/spectralwidth;

%calculate the time scale
t=[0:dwelltime:(sz(1)-1)*dwelltime];

%FILLING IN DATA STRUCTURE
out=in;
out.fids=fids;
out.specs=specs;
out.sz=sz;
out.ppm=ppm;
out.t=t;
out.spectralwidth=spectralwidth;
out.dwelltime=dwelltime;

%FILLING IN THE FLAGS
out.flags=in.flags;
out.flags.writtentostruct=1;
out.flags.freqranged=1;

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function obj = updateImage(obj, imageIndex)

% HEATMAPS
% z: ...[DONE]
% x: ...[DONE]
% y: ...[DONE]
% name: ...[DONE]
% zauto: ...[DONE]
% zmin: ...[DONE]
% zmax: ...[DONE]
% colorscale: ...[DONE]
% reversescale: ...[DONE]
% showscale: ...[DONE]
% colorbar: ...[HANDLED BY COLORBAR]
% zsmooth: ...[NOT SUPPORTED BY MATLAB]
% opacity: ---[TODO]
% xaxis: ...[DONE]
% yaxis: ...[DONE]
% showlegend: ...[DONE]
% stream: ...[HANDLED BY PLOTLYSTREAM]
% visible: ...[DONE]
% x0: ...[NOT SUPPORTED IN MATLAB]
% dx: ...[NOT SUPPORTED IN MATLAB]
% y0: ...[NOT SUPPORTED IN MATLAB]
% dy: ...[NOT SUPPORTED IN MATLAB]
% xtype: ...[NOT SUPPORTED IN MATLAB]
% ytype: ...[NOT SUPPORTED IN MATLAB]
% type: ...[DONE]

%-FIGURE STRUCTURE-%
figure_data = get(obj.State.Figure.Handle);

%-AXIS STRUCTURE-%
axis_data = get(obj.State.Plot(imageIndex).AssociatedAxis);

%-AXIS INDEX-%
axIndex = obj.getAxisIndex(obj.State.Plot(imageIndex).AssociatedAxis);

%-CHECK FOR MULTIPLE AXES-%
[xsource, ysource] = findSourceAxis(obj,axIndex);

%-IMAGE DATA STRUCTURE- %
image_data = get(obj.State.Plot(imageIndex).Handle);

%-AXIS DATA-%
eval(['xaxis = obj.layout.xaxis' num2str(xsource) ';']);
eval(['yaxis = obj.layout.yaxis' num2str(ysource) ';']);

%-------------------------------------------------------------------------%

%-image xaxis-%
obj.data{imageIndex}.xaxis = ['x' num2str(xsource)];

%-------------------------------------------------------------------------%

%-image yaxis-%
obj.data{imageIndex}.yaxis = ['y' num2str(ysource)];

%-------------------------------------------------------------------------%

%-image type-%
obj.data{imageIndex}.type = 'heatmap';

%-------------------------------------------------------------------------%

%-image x-%
if (size(image_data.XData,2) == 2)
    obj.data{imageIndex}.x = image_data.XData(1):image_data.XData(2);
else
    obj.data{imageIndex}.x = image_data.XData;
end

%-------------------------------------------------------------------------%

%-image y-%
if (size(image_data.YData,2) == 2)
    obj.data{imageIndex}.y = image_data.YData(1):image_data.YData(2);
else
    obj.data{imageIndex}.y = image_data.YData;
end

%-------------------------------------------------------------------------%

%-image z-%
if(size(image_data.CData,3) > 1)
    % TODO: ALLOW FOR TRUE COLOUR SPECS.
    obj.data{imageIndex}.z = image_data.CData(:,:,1);
else
    obj.data{imageIndex}.z = image_data.CData;
end

%-------------------------------------------------------------------------%

%-image name-%
obj.data{imageIndex}.name = image_data.DisplayName;

%-------------------------------------------------------------------------%

%-image visible-%
obj.data{imageIndex}.visible = strcmp(image_data.Visible,'on');

%-------------------------------------------------------------------------%

%-image showscale-%
obj.data{imageIndex}.showscale = false;

%-------------------------------------------------------------------------%

%-image zauto-%
obj.data{imageIndex}.zauto = false;

%-------------------------------------------------------------------------%

%-image zmin-%
obj.data{imageIndex}.zmin = axis_data.CLim(1);

%-------------------------------------------------------------------------%

%-image zmax-%
obj.data{imageIndex}.zmax = axis_data.CLim(2);

%-------------------------------------------------------------------------%

%-COLORSCALE (ASSUMES IMAGE CDATAMAP IS 'SCALED')-%

%-image colorscale-%
colormap = figure_data.Colormap;

for c = 1:length(colormap)
    col =  255*(colormap(c,:));
    obj.data{imageIndex}.colorscale{c} = {(c-1)/length(colormap), ['rgb(' num2str(col(1)) ',' num2str(col(2)) ',' num2str(col(3)) ')']};
end

%-------------------------------------------------------------------------%

%-image showlegend-%
leg = get(image_data.Annotation);
legInfo = get(leg.LegendInformation);

switch legInfo.IconDisplayStyle
    case 'on'
        showleg = true;
    case 'off'
        showleg = false;
end

obj.data{imageIndex}.showlegend = showleg;

%-------------------------------------------------------------------------%

end

author: goGPS-Project


text:
function fd1d_heat_implicit_test03 ( )

%*****************************************************************************80
%
%% FD1D_HEAT_IMPLICIT_TEST03 does a simple test problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD1D_HEAT_IMPLICIT_TEST03:\n' );
  fprintf ( 1, '  Compute an approximate solution to the time-dependent\n' );
  fprintf ( 1, '  one dimensional heat equation:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    dH/dt - K * d2H/dx2 = f(x,t)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Run a simple test case.\n' );
%
%  Heat coefficient.
%
  k = k_test03 ( );
%
%  X_NUM is the number of equally spaced nodes to use between 0 and 1.
%
  x_num = 21;
  x_min = -5.0;
  x_max = +5.0;
  dx = ( x_max - x_min ) / ( x_num - 1 );
  x = linspace ( x_min, x_max, x_num );
%
%  T_NUM is the number of equally spaced time points between 0 and 10.0.
%
  t_num = 81;
  t_min = 0.0;
  t_max = 4.0;
  dt = ( t_max - t_min ) / ( t_num - 1 );
  t = linspace ( t_min, t_max, t_num );
%
%  Get the CFL coefficient.
%
  cfl = fd1d_heat_implicit_cfl ( k, t_num, t_min, t_max, x_num, x_min, x_max );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of X nodes = %d\n', x_num );
  fprintf ( 1, '  X interval is [%f,%f]\n', x_min, x_max );
  fprintf ( 1, '  X spacing is %f\n', dx );
  fprintf ( 1, '  Number of T values = %d\n', t_num );
  fprintf ( 1, '  T interval is [%f,%f]\n', t_min, t_max );
  fprintf ( 1, '  T spacing is %f\n', dt );
  fprintf ( 1, '  Constant K = %g\n', k );
  fprintf ( 1, '  CFL coefficient = %g\n', cfl );
%
%  Get the system matrix.
%
  a = fd1d_heat_implicit_matrix ( x_num, cfl );

  hmat = zeros ( x_num, t_num );

  for j = 1 : t_num
    if ( j == 1 )
      h = ic_test03 ( x_num, x, t(j) );
      h = bc_test03 ( x_num, x, t(j), h );
    else
      h = fd1d_heat_implicit ( a, x_num, x, t(j-1), dt, cfl, @rhs_test03, @bc_test03, h );
    end
    hmat(1:x_num,j) = h(1:x_num);
  end
%
%  Plot the data.
%
  figure ( 3 )
  [ tmat, xmat ] = meshgrid ( t, x );
  mesh ( tmat, xmat, hmat );
  title ( 'H(X,T) for TEST03 computed by FD1D\_HEAT\_IMPLICIT' );
  xlabel ( '<-- Time -->' );
  ylabel ( '<-- X -->' );
  zlabel ( '<--  H(X,T) -->' );
%
%  Write the data to files.
%
  r8mat_write ( 'h_test03.txt', x_num, t_num, hmat );
  r8vec_write ( 't_test03.txt', t_num, t );
  r8vec_write ( 'x_test03.txt', x_num, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  H(X,T) written to "h_test03.txt"\n' );
  fprintf ( 1, '  T values written to "t_test03.txt"\n' );
  fprintf ( 1, '  X values written to "x_test3.txt"\n' );
  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef PTKUnclosedLungIncludingTrachea < PTKPlugin
    % PTKUnclosedLungIncludingTrachea. Plugin to segment the lung regions
    %     including the airways.
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Lungs with trachea';
        ToolTip = 'Find unclosed lung region including the trachea and main bronchi'
        Category = 'Lungs'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = true
        Visibility = 'Developer'
        Version = 3
    end

    methods (Static)
        function results = RunPlugin(dataset, reporting)
            if dataset.IsGasMRI
                results = dataset.GetResult('PTKSegmentGasMRI');
                results.AddBorder(1);
                results = PTKGetMainRegionExcludingBorder(results, 1000000, reporting);
                results.RemoveBorder(1);
            elseif strcmp(dataset.GetImageInfo.Modality, 'MR')
                lung_threshold = dataset.GetResult('PTKMRILungThreshold');
                results = lung_threshold.LungMask;

            else
                threshold_image = dataset.GetResult('PTKThresholdLungFiltered');
                threshold_image.ChangeRawImage(threshold_image.RawImage > 0);

                reporting.ShowProgress('Searching for largest connected region');

                % Find the main component, excluding any components touching the border
                original_image = dataset.GetResult('PTKLungROI');
                threshold_image = PTKGetMainRegionExcludingPaddingBorder(original_image, threshold_image, 1000000, false, reporting);

                results = threshold_image;
                results.ImageType = PTKImageType.Colormap;
            end
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [charge, chargeWithoutProtons] = getChargeFromInChI(InChI)
% Returns the charge from a given InChI string
%
% USAGE:
%
%    [charge, chargeWithoutProtons] = getChargeFromInChI(InChI)
%
% INPUT:
%    InChI:                   The Inchi Identifier - string
%
% OUTPUTS:
%    charge:                  The charge encoded in the `InChi` string (including protonation)
%    chargeWithoutProtons:    The charge encoded in the `InChi` ignoring the protonation state
%
% NOTE:
%
%    InChI Charge is defined in the charge layer and can be modified in the
%    proton layer. If nothing is defined, the compound is uncharged.
%    First: Discard any "Reconnected" parts, as those don't influence the
%    charge
%
% .. Author:
%       - Ronan Fleming, 23 Sept 09
%       - Thomas Pfau, May 2017, Updated

    InChI = regexprep(InChI,'/r.*','');

    %Charge Layer: (either at the end or at the start)
    q_layer = regexp(InChI,'/q(.*?)/|/q(.*?)$','tokens');
    %We need to store a potential change of charge in the fixed-H
    %layer, as we need to ignore that change later on.
    q_fixed_H_layer = regexp(InChI,'/f.*?/q(.*?)/|/f.*?/q(.*?)$','tokens');

    %proton layer
    p_layer = regexp(InChI,'/p(.*?)/|/p(.*?)$','tokens');

    chargeWithoutProtons = 0;

    if ~isempty(q_layer)
        %Get individual charges from splitted reactions.
        individualCharges = cellfun(@(x) {strsplit(x{1},';')},q_layer);
        %And calculate the charge by evaluating the individual components.
        chargeWithoutProtons = sum(cellfun(@(x) sum(cellfun(@(y) eval(y) , x)), individualCharges));
        if ~isempty(q_fixed_H_layer)
            individualCharges_fixed_H = cellfun(@(x) {strsplit(x{1},';')},q_fixed_H_layer);
            chargeWithoutProtons = chargeWithoutProtons - sum(cellfun(@(x) sum(cellfun(@(y) eval(y) , x)), individualCharges_fixed_H));
        end
    end

    proton_charges = 0;
    if ~isempty(p_layer)
        individualProtons = cellfun(@(x) {strsplit(x{1},';')},p_layer);
        proton_charges = cellfun(@(x) sum(cellfun(@(y) eval(y) , x)), individualProtons);
    end
    %The overall Charge is the combination of charge from protons and base
    %charge
    charge = proton_charges + chargeWithoutProtons;

author: opencobra

repo: cobratoolbox


text:
function [ labels, valid_ids, filenames  ] = extract_FERA2011_labels( FERA2011_dir, recs, aus )
%EXTRACT_SEMAINE_LABELS Summary of this function goes here
%   Detailed explanation goes here

    num_files = numel(recs);

    % speech invalidates lower face AUs
    labels = cell(num_files, 1);
    valid_ids = cell(num_files, 1);
    filenames = cell(num_files, 1);

    file_id = 1;

    for i=1:numel(recs)

        file = [FERA2011_dir, '/', recs{i}, '/', recs{i}, '-au.dat'];

        [~, filename,~] = fileparts(file);
        filenames{file_id} = filename;

        data = dlmread(file); %import annotations for one video file

        speech = data(:,end);

        labels{file_id} = data(:, aus);

        % Finding the invalid regions
        if(aus(1) >= 10)
            valid = speech == 0;
        else
            valid = true(size(speech,1), 1);
        end

        % all indices in SEMAINE are valid
        valid_ids{file_id} = valid;

        file_id = file_id + 1;
    end

    labels = labels(1:file_id-1);
    valid_ids = valid_ids(1:file_id-1);
    filenames = filenames(1:file_id-1);

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function state = isintvector(var)

if ~isempty(var) & ...
   isreal(var) & ...
   isequal(var,fix(var)) & ...
   ndims(var) == 2 & ...
   ~(size(var,1)>1 & size(var,2)>1)
   state = 1;
else
   state = 0;
end

repo: mcbench-benchmarks


text:
function SNR = hsvargplvmShowSNR(model, layers, displ)

if nargin < 3
    displ = true;
end

if nargin < 2 || isempty(layers)
    layers = 1:model.H;
end

for h=layers
    if displ
        fprintf('# SNR Layer %d\n',h)
    end
    for m=1:model.layer{h}.M
        if isfield(model.layer{h}.comp{m}, 'mOrig')
            varY = var(model.layer{h}.comp{m}.mOrig(:));
        else
            varY = var(model.layer{h}.comp{m}.m(:));
        end
        beta = model.layer{h}.comp{m}.beta;
        SNR{h}(m) = varY * beta;
        if displ
            fprintf('    Model %d: %f  (varY=%f, 1/beta=%f)\n', m, SNR{h}(m), varY, 1/beta)
        end
    end
end

author: SheffieldML

repo: deepGP

save_path: github-repos/MATLAB/SheffieldML-deepGP


text:
function create_table_of_algorithms_comb(up)
%CREATE_TABLE_OF_ALGORITHMS creates a table of algorithms and their performances
%	            create_table_of_algorithms(up)
%
%	Inputs:
%		data            data files should be stored in the specified format
%       up              universal parameters structure
%
%	Outputs:
%       tables          csv tables of results
%

fprintf('\n--- Creating Results Tables ');

% save setup
save_name = up.paths.filenames.results_table;
savepath = [up.paths.tables_save_folder, save_name, '.mat'];
if ~up.analysis.redo_stats
    exist_log = check_exists(savepath, save_name);
    if exist_log
        return
    end
end

% Load BA Data
load_name = 'BA_results';
loadpath = [up.paths.data_save_folder, up.paths.filenames.global_BA, '.mat'];
load(loadpath, load_name);

% Load algorithm names
load_name = up.paths.filenames.alg_names;
loadpath = [up.paths.data_save_folder, up.paths.filenames.alg_names, '.mat'];
load(loadpath, load_name);

% Load study data
load_name = 'trad_stats';
loadpath = [up.paths.data_save_folder, up.paths.filenames.study_stats, '.mat'];
load(loadpath, load_name);

% identify groups
groups = fieldnames(BA_results);

%% Create a BA table and a Stats table for each group

for group_no = 1 : length(groups)
    group = groups{group_no};

    %% Extract relevant data
    rel_data = extract_rel_data(group, BA_results, trad_stats, up);

    %% Make Table of BA Results
    BA_results_table = create_table_results(alg_names, rel_data, up, group, 'BA');

    %% Save to file
    save_name = [up.paths.filenames.BA_results_table, '_', group];
    savepath = [up.paths.tables_save_folder, save_name];
    xlswrite(savepath, BA_results_table)

    %% Make Table of Stats Results
    stats_results_table = create_table_results(alg_names, rel_data, up, group, 'stats');

    %% Save to file
    save_name = [up.paths.filenames.stats_results_table, '_', group];
    savepath = [up.paths.tables_save_folder, save_name];
    xlswrite(savepath, stats_results_table)

    %% Make combined table of results
    comb_results_table = create_table_results(alg_names, rel_data, up, group, 'comb');

    %% Save to file
    save_name = [up.paths.filenames.comb_results_table, '_', group];
    savepath = [up.paths.tables_save_folder, save_name];
    xlswrite(savepath, comb_results_table)

end

end

function rel_data = extract_rel_data(group, BA_results, stats, up)

rel_data.BA = extractfield(BA_results, group); rel_data.BA = rel_data.BA{1,1};
rel_data.stats = extractfield(stats, group); rel_data.stats = rel_data.stats{1,1};

end

function table_data = create_table_results(alg_names, rel_data, up, group, type)

% Find alg numbers:
alg_names.no = cell(length(alg_names.names),1);
for s = 1 : length(alg_names.no)
    alg_names.no{s} = ['a', num2str(s)];
end

% Extract relevant variables for table:
tab_vars.alg_no = alg_names.no;
if isfield(alg_names.meths, 'xa')
    tab_vars.m_xa = alg_names.meths.xa;
    tab_vars.m_xb = alg_names.meths.xb;
    tab_vars.m_ef = alg_names.meths.ef;
    tab_vars.m_et = alg_names.meths.et;
    tab_vars.m_fm = alg_names.meths.fm;
    tab_vars.m_ft = alg_names.meths.ft;
    % modify the fme:
    rel_els = ~isnan(tab_vars.m_fm);
    tab_vars.m_xb(rel_els) = 99;
end
tab_vars.alg_name = alg_names.names;
tab_vars.signal = alg_names.sigs;

% Extract relevant statistics for table:
if strcmp(type, 'BA')
    tab_vars.bias = rel_data.BA.bias.val;
    tab_vars.two_sd = rel_data.BA.two_sd.val;
    tab_vars.prop_wins_est = rel_data.stats.prop_wins_est;
    tab_vars.prop_wins_ref_and_good_sqi = rel_data.stats.prop_wins_good_sqi_and_ref;
elseif strcmp(type, 'stats')
    tab_vars.percerr = rel_data.stats.percerr;
    tab_vars.mae = rel_data.stats.mae;
    tab_vars.sdae = rel_data.stats.sdae;
    tab_vars.rmse = rel_data.stats.rmse;
    tab_vars.cp2 = rel_data.stats.cp2;
    tab_vars.prop_wins_est = rel_data.stats.prop_wins_est;
    tab_vars.total_wins_all = rel_data.stats.total_wins_all;
elseif strcmp(type, 'comb')
    tab_vars.two_sd = rel_data.BA.two_sd.val;
    tab_vars.bias = rel_data.BA.bias.val;
    tab_vars.cp2 = rel_data.stats.cp2;
    tab_vars.percerr = rel_data.stats.percerr;
    tab_vars.mae = rel_data.stats.mae;
    tab_vars.sdae = rel_data.stats.sdae;
    tab_vars.rmse = rel_data.stats.rmse;
    tab_vars.prop_wins_ref_and_good_sqi = rel_data.stats.prop_wins_good_sqi_and_ref;
    tab_vars.prop_wins_est = rel_data.stats.prop_wins_est;
    tab_vars.total_wins_all = rel_data.stats.total_wins_all;
end

%% Determine what goes into the table
headers = fieldnames(tab_vars);
table_data = cell(1+length(tab_vars.alg_name), length(headers));
table_data(1,:) = headers;
for col_no = 1 : length(headers)
    eval(['col_data = tab_vars.' headers{col_no} ';']);
    if ~iscell(col_data)
        table_data(2:end,col_no) = num2cell(col_data);
    else
        table_data(2:end,col_no) = col_data;
    end
end

end

function [subj_data, BA_data, stats_data, old_names] = load_data(group, up)

% Load subj data
load_name = ['data_', group];
loadpath = [up.paths.data_save_folder, up.paths.filenames.group_data, group, '.mat'];
load(loadpath, load_name)
eval(['subj_data = data_' group  ';']);

% Load B-A results
load_name = 'BA';
loadpath = [up.paths.data_save_folder, up.paths.filenames.group_BA, group, '.mat'];
load(loadpath, load_name)
BA_data = BA; clear BA

% Load Stats data
load_name = 'acc';
loadpath = [up.paths.data_save_folder, up.paths.filenames.group_stats, group, '.mat'];
load(loadpath, load_name)
stats_data = acc;
old_names = acc.rr_alg_names;

end

function create_table_stats_results_subj(alg_names, sig_names, old_names, subj_data, up, group)

%% Find prop with an abs error of < 1 bpm
eval(['subjs = unique(subj_data.' old_names{1} '.subj);']);
prop_acc = zeros(length(old_names), length(subjs));
for alg_no = 1 : length(old_names)
     eval(['rel_data = subj_data.' old_names{alg_no} ';']);

     rel_data.error = rel_data.est - rel_data.ref;
     rel_data.abs_error = abs(rel_data.error);
     rel_data.acc_log = rel_data.abs_error < 1;

     for subj_no = 1:length(subjs)
         subj = subjs(subj_no);
         rel_els = logical(rel_data.subj == subj);
         prop_acc(alg_no, subj_no) = sum(rel_data.acc_log(rel_els)) / length(rel_data.acc_log(rel_els));
     end

     % extract relevant parts of algs:
     for comp_str = {'FMe', 'Flt', 'RSa', 'PDt', 'Est', 'MFu', 'TFu'}
         if ~exist([lower(comp_str{1,1}), '_meth'], 'var')
             eval([lower(comp_str{1,1}), '_meth = nan(length(old_names),1);']);
         end
         rel_el = strfind(alg_names{alg_no}, comp_str{1,1});
         if rel_el
             temp = str2num(alg_names{alg_no}(rel_el+3));
             eval([lower(comp_str{1,1}), '_meth(alg_no) = temp;']);
         end
     end

     clear rel_data
end

%% Determine what goes into the table
headers = {'Signal', 'Algorithm Name', 'Old Name', 'FMe', 'Flt', 'RSa', 'PDt', 'Est', 'MFu', 'TFu'};
for subj_no = 1:length(subjs)
    headers{end+1} = ['Prop_acc ' num2str(subjs(subj_no))];
end

table_data = cell(1+length(alg_names), length(headers));
table_data(1,:) = headers;
table_data(2:end,1) = sig_names;
table_data(2:end,2) = alg_names;
table_data(2:end,3) = old_names;
table_data(2:end,4) = num2cell(fme_meth);
table_data(2:end,5) = num2cell(flt_meth);
table_data(2:end,6) = num2cell(rsa_meth);
table_data(2:end,7) = num2cell(pdt_meth);
table_data(2:end,8) = num2cell(est_meth);
table_data(2:end,9) = num2cell(mfu_meth);
table_data(2:end,10) = num2cell(tfu_meth);
for s = 1 : size(prop_acc,2)
    table_data(2:end,10+s) = num2cell(prop_acc(:,s));
end

%% Save to file
save_name = [up.paths.filenames.subj_results_table '_' group];
savepath = [up.paths.tables_save_folder, save_name];
xlswrite(savepath,table_data)

end

function create_table_BA_results(alg_names, sig_names, BA_data, up, group)

%% Determine what goes into the table
headers = {'Signal', 'Algorithm Name', 'Bias', 'Bias (LCI)', 'Bias (UCI)', 'LLOA', 'LLOA (LCI)', 'LLOA (UCI)', 'ULOA', 'ULOA (LCI)', 'ULOA (UCI)', 'Precision', 'Precision (LCI)', 'Precision (UCI)'};
table_data = cell(1+length(alg_names), length(headers));
table_data(1,:) = headers;
table_data(2:end,1) = sig_names;
table_data(2:end,2) = alg_names;
table_data(2:end,3) = num2cell(BA_data.bias.val);
table_data(2:end,4) = num2cell(BA_data.bias.lci);
table_data(2:end,5) = num2cell(BA_data.bias.uci);
table_data(2:end,6) = num2cell(BA_data.lloa.val);
table_data(2:end,7) = num2cell(BA_data.lloa.lci);
table_data(2:end,8) = num2cell(BA_data.lloa.uci);
table_data(2:end,9) = num2cell(BA_data.uloa.val);
table_data(2:end,10) = num2cell(BA_data.uloa.lci);
table_data(2:end,11) = num2cell(BA_data.uloa.uci);
table_data(2:end,12) = num2cell(BA_data.prec.val);
table_data(2:end,13) = num2cell(BA_data.prec.lci);
table_data(2:end,14) = num2cell(BA_data.prec.uci);

%% Save to file
save_name = [up.paths.filenames.BA_results_table '_' group];
savepath = [up.paths.tables_save_folder, save_name];
xlswrite(savepath,table_data)

end

function create_table_stats_results(alg_names, sig_names, stats_data, up, group)

%% Determine what goes into the table
headers = {'Signal', 'Algorithm Name', 'MAE', 'MAE (lq)', 'MAE (uq)', 'SDAE', 'SDAE (lq)', 'SDAE (uq)', 'PERC ERR', 'PERC ERR (lq)', 'PERC ERR (uq)', 'RMSE', 'RMSE (lq)', 'RMSE (uq)', 'MEAN ERROR', 'MEAN ERROR (lq)', 'MEAN ERROR (uq)', 'STD ERROR', 'STD ERROR (lq)', 'STD ERROR (uq)', 'PROP WINS', 'PROP WINS (lq)', 'PROP WINS (uq)'};
table_data = cell(1+length(alg_names), length(headers));
table_data(1,:) = headers;
table_data(2:end,1) = sig_names;
table_data(2:end,2) = alg_names;
table_data(2:end,3) = num2cell(stats_data.mae.med);
table_data(2:end,4) = num2cell(stats_data.mae.lq);
table_data(2:end,5) = num2cell(stats_data.mae.uq);
table_data(2:end,6) = num2cell(stats_data.sdae.med);
table_data(2:end,7) = num2cell(stats_data.sdae.lq);
table_data(2:end,8) = num2cell(stats_data.sdae.uq);
table_data(2:end,9) = num2cell(stats_data.percerr.med);
table_data(2:end,10) = num2cell(stats_data.percerr.lq);
table_data(2:end,11) = num2cell(stats_data.percerr.uq);
table_data(2:end,12) = num2cell(stats_data.rmse.med);
table_data(2:end,13) = num2cell(stats_data.rmse.lq);
table_data(2:end,14) = num2cell(stats_data.rmse.uq);
table_data(2:end,15) = num2cell(stats_data.mean_diff.med);
table_data(2:end,16) = num2cell(stats_data.mean_diff.lq);
table_data(2:end,17) = num2cell(stats_data.mean_diff.uq);
table_data(2:end,18) = num2cell(stats_data.std_diff.med);
table_data(2:end,19) = num2cell(stats_data.std_diff.lq);
table_data(2:end,20) = num2cell(stats_data.std_diff.uq);
table_data(2:end,21) = num2cell(stats_data.prop_wins.med);
table_data(2:end,22) = num2cell(stats_data.prop_wins.lq);
table_data(2:end,23) = num2cell(stats_data.prop_wins.uq);

%% Save to file
save_name = [up.paths.filenames.stats_results_table '_' group];
savepath = [up.paths.tables_save_folder, save_name];
xlswrite(savepath,table_data)

end

author: peterhcharlton

repo: RRest


text:
function [A, b, x0, lambda, lambda_max] = generate_lasso_data(n, d, k, noise_level)

    [A,~] = qr(randn(n,d),0);
    A = A';
    p = randperm(n);
    p = p(1:k);                                 % select location of k nonzeros
    x0 = zeros(n,1);
    x0(p) = randn(k,1);
    b = A*x0 + noise_level*randn(d, 1);                  % add random noise
    lambda_max = norm( A'*b, 'inf' );
    lambda = 0.1*lambda_max;

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function asa066_test02 ( )

%*****************************************************************************80
%
%% TEST02 tests NORMP.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02:\n' );
  fprintf ( 1, '  Compare tabulated values of the normal\n' );
  fprintf ( 1, '  Cumulative Density Function against values\n' );
  fprintf ( 1, '  computed by NORMP.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '         X        CDF                       CDF' );
  fprintf ( 1, '                    DIFF\n' );
  fprintf ( 1, '               (tabulated)                 (NORMP)\n' );
  fprintf ( 1, '\n' );

  upper = 0;
  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = normal_01_cdf_values ( n_data );

    if ( n_data == 0 )
      break
    end

    [ p, q, pdf ] = normp ( x );
    fx2 = p;

    fprintf ( 1, '  %10.4e  %24.16e  %24.16e  %10.4e\n', ...
    x, fx, fx2, abs ( fx - fx2 ) );

  end

  return
end


text:
function PlotAggregationVariance(AggregationVariance,SimulSecMoms,SecMom_Empirical,...
    SecMom_Annualized,NumSimulations,b,s_2,Name)
figure
for n=1:NumSimulations
    hold on
    h=plot(AggregationVariance,SimulSecMoms(:,n),'.');
    set(h,'color','k','markersize',1)
end
hold on
h=plot(AggregationVariance,SecMom_Empirical);
set(h,'color','r','linewidth',2);
hold on
h=plot(AggregationVariance,SecMom_Annualized);
set(h,'color','k','linewidth',2);
set(gca,'ytick',[],'xlim',[AggregationVariance(1) AggregationVariance(end)])
grid on
xlabel('aggregation size (days)')
ylabel('aggregated variance')
set(gcf,'Name',['               long memory of ' Name ' par swap rate'])

figure
h2=plot(log(AggregationVariance),log(SecMom_Empirical),'.');
%set(h,'color','k');
hold on
%h2=plot(log(AggregationVariance),b(1)+b(2)*log(AggregationVariance),'r');
hold on
y_norm=log(s_2)+log(AggregationVariance);
h3=plot(log(AggregationVariance),y_norm,'k');
grid on
legend([h2 h3],['frac. B. m.: H=' num2str(b(2)/2)],'B. m.: H=0.5','location','northwest')
xlabel('log-aggregation size (days)')
ylabel('log-aggregated variance')
set(gcf,'Name',['               long memory of ' Name ' par swap rate'])


text:
function varargout = configure_classify(varargin)
% CONFIGURE_CLASSIFY M-file for configure_classify.fig
%      CONFIGURE_CLASSIFY, by itself, creates a new CONFIGURE_CLASSIFY or raises the existing
%      singleton*.
%
%      H = CONFIGURE_CLASSIFY returns the handle to a new CONFIGURE_CLASSIFY or the handle to
%      the existing singleton*.
%
%      CONFIGURE_CLASSIFY('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in CONFIGURE_CLASSIFY.M with the given input arguments.
%
%      CONFIGURE_CLASSIFY('Property','Value',...) creates a new CONFIGURE_CLASSIFY or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before configure_classify_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to configure_classify_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

% Copyright 2002-2003 The MathWorks, Inc.

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

% Last Modified by GUIDE v2.5 01-May-2006 23:22:44

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @configure_classify_OpeningFcn, ...
                   'gui_OutputFcn',  @configure_classify_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 configure_classify is made visible.
function configure_classify_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 configure_classify (see VARARGIN)

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

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes configure_classify wait for user response (see UIRESUME)
% uiwait(handles.figure1);
% configure_classify(0,7500,400,500,44100,[0.01 0.001],[3 5],[0 20000])

handles.lowerfreq = varargin{1};  % Lower frequency for zooming
handles.upperfreq = varargin{2}; % Upper frequency for zooming

handles.classified_height = varargin{3} ; % the height of the image in the classified axes
handles.classified_width = varargin{4};  % the width of the image in the classified axes

handles.Fs = varargin{5}; % Frequency of audio sampling per second
handles.movingwin = varargin{6}; % Size of the moving window in seconds; the first number is the window size and the second is the step size
handles.tapers = varargin{7}; % Tapers for smoothing
handles.fpass = varargin{8}; % Range of frequency sampling

handles.fixed = varargin{9}; % Fixed scaling of the classified axes

set(handles.ZoomLowerFreq,'String',num2str(handles.lowerfreq));
set(handles.ZoomUpperFreq,'String',num2str(handles.upperfreq));
set(handles.ClassifiedWidth,'String',num2str(handles.classified_width));
set(handles.ClassifiedHeight,'String',num2str(handles.classified_height));
set(handles.Frequency,'String',num2str(handles.Fs));
set(handles.WinSize,'String',num2str(handles.movingwin(1) * 1000));
set(handles.StepSize,'String',num2str(handles.movingwin(2) * 1000));
set(handles.TW,'String',num2str(handles.tapers(1)));
set(handles.MinFreq,'String',num2str(handles.fpass(1)));
set(handles.MaxFreq,'String',num2str(handles.fpass(2)));

set(handles.FixedCheckbox,'Value',handles.fixed);

% set(handles.ZoomLowerFreq,'Enable','off');
% set(handles.ZoomUpperFreq,'Enable','off');
% set(handles.ClassifiedWidth,'Enable','off');
% set(handles.ClassifiedHeight,'Enable','off');
% set(handles.Frequency,'Enable','off');
% set(handles.WinSize,'Enable','off');
% set(handles.StepSize,'Enable','off');
% set(handles.TW,'Enable','off');
% set(handles.MinFreq,'Enable','off');
% set(handles.MaxFreq,'Enable','off');

uiwait(handles.figure1);

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

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

close;

function WinSize_Callback(hObject, eventdata, handles)
% hObject    handle to WinSize (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 WinSize as text
%        str2double(get(hObject,'String')) returns contents of WinSize as a double

% --- Executes during object creation, after setting all properties.
function WinSize_CreateFcn(hObject, eventdata, handles)
% hObject    handle to WinSize (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
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

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

lowerfreq = str2num(get(handles.ZoomLowerFreq,'String'));  % Lower frequency for zooming
upperfreq = str2num(get(handles.ZoomUpperFreq,'String')); % Upper frequency for zooming

classified_height = str2num(get(handles.ClassifiedHeight,'String')); % the height of the image in the classified axes
classified_width = str2num(get(handles.ClassifiedWidth,'String'));  % the width of the image in the classified axes

Fs = str2num(get(handles.Frequency,'String')); % Frequency of audio sampling per second

winsizeS = str2num(get(handles.WinSize,'String')) / 1000;
stepS = str2num(get(handles.StepSize,'String')) / 1000;
movingwin = [winsizeS stepS]; % Size of the moving window in seconds; the first number is the window size and the second is the step size

tw = str2num(get(handles.TW,'String'));

fpasslower = str2num(get(handles.MinFreq,'String')); % Range of frequency sampling
fpassupper = str2num(get(handles.MaxFreq,'String'));
fpass = [fpasslower fpassupper];

ierror = 1; % Indicates no errors encountered

if isempty(classified_height) || (classified_height < 1)
    ierror = 0;
end

if isempty(tw) || tw < 0
    ierror = 0;
end

if isempty(lowerfreq) || lowerfreq < 0
    ierror = 0;
end

if isempty(fpasslower) || fpasslower < 0
    ierror = 0;
end

if isempty(fpassupper) || fpassupper < fpasslower
    ierror = 0;
end

if isempty(upperfreq) || lowerfreq > upperfreq
    ierror = 0;
end

if isempty(winsizeS) || winsizeS < 0
    ierror = 0;
end

if isempty(stepS) || stepS < 0
    ierror = 0
end

if isempty(tw) || tw < 0
    ierror = 0
else
    tapers = [tw,floor(2*tw-1)]; % Tapers for smoothing
end

fixed = get(handles.FixedCheckbox,'Value');

if ierror == 0
    ;
else
    handles.output = {lowerfreq,upperfreq,classified_height,classified_width,Fs,movingwin,tapers,fpass,fixed};
    guidata(hObject,handles);
    uiresume(handles.figure1);
end

%uiresume;
%close;

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

handles.output = 0;
guidata(hObject,handles);
uiresume(handles.figure1);

function StepSize_Callback(hObject, eventdata, handles)
% hObject    handle to StepSize (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 StepSize as text
%        str2double(get(hObject,'String')) returns contents of StepSize as a double

% --- Executes during object creation, after setting all properties.
function StepSize_CreateFcn(hObject, eventdata, handles)
% hObject    handle to StepSize (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
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function TW_Callback(hObject, eventdata, handles)
% hObject    handle to TW (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 TW as text
%        str2double(get(hObject,'String')) returns contents of TW as a double

% --- Executes during object creation, after setting all properties.
function TW_CreateFcn(hObject, eventdata, handles)
% hObject    handle to TW (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
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function MinFreq_Callback(hObject, eventdata, handles)
% hObject    handle to MinFreq (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 MinFreq as text
%        str2double(get(hObject,'String')) returns contents of MinFreq as a double

% --- Executes during object creation, after setting all properties.
function MinFreq_CreateFcn(hObject, eventdata, handles)
% hObject    handle to MinFreq (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
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function MaxFreq_Callback(hObject, eventdata, handles)
% hObject    handle to MaxFreq (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 MaxFreq as text
%        str2double(get(hObject,'String')) returns contents of MaxFreq as a double

% --- Executes during object creation, after setting all properties.
function MaxFreq_CreateFcn(hObject, eventdata, handles)
% hObject    handle to MaxFreq (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
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function ZoomLowerFreq_Callback(hObject, eventdata, handles)
% hObject    handle to ZoomLowerFreq (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 ZoomLowerFreq as text
%        str2double(get(hObject,'String')) returns contents of ZoomLowerFreq as a double

% --- Executes during object creation, after setting all properties.
function ZoomLowerFreq_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ZoomLowerFreq (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
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function ZoomUpperFreq_Callback(hObject, eventdata, handles)
% hObject    handle to ZoomUpperFreq (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 ZoomUpperFreq as text
%        str2double(get(hObject,'String')) returns contents of ZoomUpperFreq as a double

% --- Executes during object creation, after setting all properties.
function ZoomUpperFreq_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ZoomUpperFreq (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
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function ClassifiedWidth_Callback(hObject, eventdata, handles)
% hObject    handle to ClassifiedWidth (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 ClassifiedWidth as text
%        str2double(get(hObject,'String')) returns contents of ClassifiedWidth as a double

% --- Executes during object creation, after setting all properties.
function ClassifiedWidth_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ClassifiedWidth (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
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function ClassifiedHeight_Callback(hObject, eventdata, handles)
% hObject    handle to ClassifiedHeight (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 ClassifiedHeight as text
%        str2double(get(hObject,'String')) returns contents of ClassifiedHeight as a double

% --- Executes during object creation, after setting all properties.
function ClassifiedHeight_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ClassifiedHeight (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
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

function Frequency_Callback(hObject, eventdata, handles)
% hObject    handle to Frequency (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 Frequency as text
%        str2double(get(hObject,'String')) returns contents of Frequency as a double

% --- Executes during object creation, after setting all properties.
function Frequency_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Frequency (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
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on button press in FixedCheckbox.
function FixedCheckbox_Callback(hObject, eventdata, handles)
% hObject    handle to FixedCheckbox (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 FixedCheckbox

handles.fixed = get(hObject,'Value');
guidata(gcbo,handles);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function out = mergeScore(x)

s1 = var(x(x>mean(x)));
s2 = var(x(x<mean(x)));

mu1 = mean(x(x>mean(x)));
mu2 = mean(x(x<mean(x)));
p  = mean(x>mean(x));

logp = zeros(numel(x), 2);

rs = zeros(numel(x), 2);
rs(x<0, 1) = 1;
rs(x>0, 2) = 1;

for k = 1:20
    logp(:,1) = -1/2*log(s1) - (x-mu1).^2/(2*s1) + log(p);
    logp(:,2) = -1/2*log(s2) - (x-mu2).^2/(2*s2) + log(1-p);

    lMax = max(logp,[],2);
    logp = logp - lMax;

    rs = exp(logp);

    pval = log(sum(rs,2)) + lMax;
    logP(k) = mean(pval);

    rs = rs./sum(rs,2);

    p = mean(rs(:,1));
    mu1 = (rs(:,1)' * x )/sum(rs(:,1));
    mu2 = (rs(:,2)' * x )/sum(rs(:,2));

    s1 = (rs(:,1)' * (x-mu1).^2 )/sum(rs(:,1));
    s2 = (rs(:,2)' * (x-mu2).^2 )/sum(rs(:,2));

end

ilow = rs(:,1)>rs(:,2);

plow = mean(rs(ilow,1));
phigh = mean(rs(~ilow,2));

% when do I split
out =  ~(plow>.9 && phigh>.9);

% if sign(mu1*mu2)>0
%     out = 0;
% end

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function save_joint_pair_stats( expidx )

p = exp_params(expidx);

pwIdxsAllrel = build_joint_pairs(p);

num_relations = length(pwIdxsAllrel);

graph = zeros(num_relations, 2);

for k = 1:num_relations
    graph(k, :) = pwIdxsAllrel{k};
end

graph = [graph; graph(:,2), graph(:,1)];

save_dir = fullfile(p.exp_dir, 'data');
load(fullfile(save_dir, 'all_pairs_stats'), 'means', 'std_devs');

means = [means; -means];
std_devs = [std_devs; std_devs];

input_file = fullfile(save_dir, 'all_stats.txt');

fid = fopen(input_file, 'wt');
write_text_matrix(fid, graph, 'graph');
write_text_matrix(fid, means, 'means');
write_text_matrix(fid, std_devs, 'std_devs');
fclose(fid);

save(fullfile(save_dir, 'all_pairs_stats_all'), 'means', 'std_devs', 'graph');

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [ o, x, w ] = cn_leg_01_1 ( n )

%*****************************************************************************80
%
%% CN_LEG_01_1 implements the midpoint rule for region CN_LEG.
%
%  Discussion:
%
%    The rule has order O = 1.
%
%    The rule has precision P = 1.
%
%    CN_LEG is the cube [-1,+1]^N with the Legendre weight function
%
%      w(x) = 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%
%    Input, integer O, the order.
%
%    Output, real X(N,O), the abscissas.
%
%    Output, real W(O), the weights.
%
  o = 1;

  expon = 0;
  value1 = c1_leg_monomial_integral ( expon );
  volume = value1 ^ n;

  x = zeros ( n, o );
  w = zeros ( o, 1 );

  k = 0;
%
%  1 point.
%
  k = k + 1;
  w(k) = volume;

  return
end


text:
close all;
clear all;
clc;
png_export = true;
pdf_export = false;

load('bin/rand_dict_bp_omp_success_with_k_figure_1.mat');

mf = spx.graphics.Figures();

mf.new_figure('Recovery probability with K');
hold all;
plot(Ks, bp_success_with_k);
plot(Ks, omp_success_with_k);
xlabel('Sparsity level');
ylabel('probability of successful recovery');
legend('BP', 'OMP');
grid on;

author: indigits

repo: sparse-plex


text:
function [ ival, pint ] = plane_normal_line_exp_int_3d ( pp, normal, p1, p2 )

%*****************************************************************************80
%
%% PLANE_NORMAL_LINE_EXP_INT_3D: intersection of plane and line in 3D.
%
%  Discussion:
%
%    The normal form of a plane in 3D is:
%
%      PP is a point on the plane,
%      N is a normal vector to the plane.
%
%    The explicit form of a line in 3D is:
%
%      P1, P2 are two points on the line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real PP(3), a point on the plane.
%
%    Input, real NORMAL(3), a normal vector to the plane.
%
%    Input, real P1(3), P2(3), two distinct points on the line.
%
%    Output, integer IVAL, the kind of intersection;
%    0, the line and plane seem to be parallel and separate;
%    1, the line and plane intersect at a single point;
%    2, the line and plane seem to be parallel and joined.
%
%    Output, real PINT(3), the coordinates of a
%    common point of the plane and line, when IVAL is 1 or 2.
%
  dim_num = 3;
%
%  Make sure the line is not degenerate.
%
  if ( line_exp_is_degenerate_nd ( dim_num, p1, p2 ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PLANE_NORMAL_LINE_EXP_INT_3D - Fatal error!\n' );
    fprintf ( 1, '  The line is degenerate.\n' );
    error ( 'PLANE_NORMAL_LINE_EXP_INT_3D - Fatal error!' );
  end
%
%  Make sure the plane normal vector is a unit vector.
%
  temp = sqrt ( sum ( normal(1:dim_num).^2 ) );

  if ( temp == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PLANE_NORMAL_LINE_EXP_INT_3D - Fatal error!\n' );
    fprintf ( 1, '  The normal vector of the plane is degenerate.\n' );
    error ( 'PLANE_NORMAL_LINE_EXP_INT_3D - Fatal error!' );
  end

  normal(1:dim_num) = normal(1:dim_num) / temp;
%
%  Determine the unit direction vector of the line.
%
  direction(1:dim_num) = p2(1:dim_num) - p1(1:dim_num);
  temp = sqrt ( sum ( direction(1:dim_num).^2 ) );
  direction(1:dim_num) = direction(1:dim_num) / temp;
%
%  If the normal and direction vectors are orthogonal, then
%  we have a special case to deal with.
%
  if ( normal(1:dim_num) * direction(1:dim_num)' == 0.0 )

    temp = normal(1:dim_num) * ( p1(1:dim_num) - pp(1:dim_num) )';

    if ( temp == 0.0 )
      ival = 2;
      pint(1:dim_num) = p1(1:dim_num);
    else
      ival = 0;
      pint(1:dim_num) = Inf;
    end

    return
  end
%
%  Determine the distance along the direction vector to the intersection point.
%
  temp = ( ( pp(1:dim_num) - p1(1:dim_num) ) * normal(1:dim_num)' )...
    /  ( direction(1:dim_num) * normal(1:dim_num)' );

  ival = 1;
  pint(1:dim_num) = p1(1:dim_num) + temp * direction(1:dim_num);

  return
end


text:
%%*****************************************************************
%% randmat: generate an mxn matrix using matlab's
%%          rand or randn functions using state = k.
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

function v = randmat(m,n,k,randtype)

try
   s = rng;
   rng(k);
   if strcmp(randtype,'n')
      v = randn(m,n);
   elseif strcmp(randtype,'u')
      v = rand(m,n);
   end
   rng(s);
catch
   if strcmp(randtype,'n')
      s = randn('state');
      randn('state',k);
      v = randn(m,n);
      randn('state',s);
   elseif strcmp(randtype,'u')
      s = rand('state');
      rand('state',k);
      v = rand(m,n);
      rand('state',s);
   end
end
%%******************************************************

author: yu-jiang

repo: radpbook

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


text:
function out = binsumImage(in, bin)
% Rebins image by summing bin(1)xbin(2) pixels together along x and y
% dimension, respectively. Starts at upper left corner. Rebinned image will
% have floor(size(in)./bin) pixels.
%
% out = binsum(in, bin)
di=0;
if strcmp(class(in),'dip_image')
    di=1;
    in=double(in);
end

npix=ceil(size(in)./bin);
out=zeros(npix);
for xi=1:npix(1)
    startx=(xi-1)*bin(1)+1;
    endx=min(startx+bin(1)-1,size(in,1));
    for yi=1:npix(2)
        starty=(yi-1)*bin(2)+1;
        endy=min(starty+bin(2)-1,size(in,2));
        binsum=sum(sum(in(startx:endx,starty:endy)));
        out(xi,yi)=binsum;
    end
end

if di
    out=dip_image(out);
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function mrvBrowseSVN(fileName)
%
% mrvBrowseSVN(fileName)
%
% Deprecated function - wraps vistaBrowseGit.m
%
% See also:
%   vistaBrowseGit.m
%
% (C) Stanford University, 2015
%

if nargin == 0
    help(mfilename)
    return
end

disp('This function is deprecated. Using vistaBrowseGit instead.');
vistaBrowseGit(fileName)

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function diagnostic = bisection_core(Constraints,Objective,options)
%BISECTION Core engine

% Note, outer layer added to enable arbitrary sign on ojective.
% Code here initially assumed maximization of scalar, so now we tweak the
% code a bit with swtchedsign to run any case. Fugly.

diagnostic.yalmiptime = 0;
diagnostic.solvertime = 0;
diagnostic.info = '';
diagnostic.problem = 0;

the_sign = 1;
if options.bisection.switchedsign
    the_sign = -1;
end

if length(getvariables(Objective)) > 1
    diagnostic.problem = -4;
    return
end

if  ~(isequal(getbase(Objective),[0 1]))
    diagnostic.problem = -4;
    return
end

if nargin < 3 || isempty(options)
    options = sdpsettings;
end

% Initialize the lower bound
% Typically a good guess
lower = 0;

% Silly upper bound
bestUpper = inf;

% Create an optimizer object which solves feasibility problem for a
% particular value of the sought variable
x = recover(setdiff(depends(Constraints),depends(Objective)));
if isempty(options) || isequal(options.solver,'')
    % We must figure out suitable solver
    [~,~,~,model] = export(replace(Constraints,Objective,pi));
    options.solver = model.solver.tag;
end
P = optimizer(Constraints,[],options,Objective,x);

% Safe-guard aganst using LMILAB
if ~isempty(strfind(struct(P).model.options.solver,'lmilab'))
    disp('Selected solver LMILAB lacks features making it unsuitable for BISECTION')
    disp('See why here https://yalmip.github.io/solver/lmilab/');
    disp('Select another SDP solver or/and install a better SDP solver');
    error('Cannot proceed due to poor SDP solver');
    return
end

if options.verbose;
    disp(['Selected solver: ' options.solver]);
    fprintf(['Generating initial bound: ' num2str(lower*the_sign)]);
end

% Make sure we actually can solve the lower problem
solvertime = tic;
[sol, flag] = P{lower};
working_sol = [];
diagnostic.solvertime = diagnostic.solvertime + toc(solvertime);
if flag == 1
    % This was infeasible, hernce we can use it as an upper bound
    i = 1;
    while flag
        bestUpper = lower;
        lower = lower - 2^(-4+i);i = i+1;
        if options.verbose;
            fprintf([' (fail), ' num2str(lower*the_sign)]);
        end
        try
            solvertime = tic;
            [sol, flag] = P{lower};
            diagnostic.solvertime = diagnostic.solvertime + toc(solvertime);
            if lower < -1e6
                  if options.verbose;
                     fprintf([' (fail). Giving up, never feasible\n']);
                  end
                diagnostic.problem = 21;
                diagnostic.info = yalmiperror(diagnostic.problem,'BISECTION');
                return
            end
        catch
            diagnostic.problem = 1;
            diagnostic.info = yalmiperror(diagnostic.problem,'BISECTION');
            return
        end
    end
    if options.verbose;
        fprintf([' (ok).']);
    end
elseif flag == 0
   working_sol = sol;
else
	diagnostic.problem = flag;
	diagnostic.info = yalmiperror(diagnostic.problem,'BISECTION');
    return
end

v = sol;
optimal = lower;
upper = bestUpper;

if isinf(upper)
    upper = lower+1;
    if options.verbose;
       fprintf([' (ok), ' num2str(upper*the_sign)]);
    end
    [sol, flag] = P{upper};
    i = 1;
    while ~flag
        % This was feasible, so we can use it as new lower bound
        lower = upper;working_sol = sol;optimal = upper;
        % Increase upper bound
        upper = upper + 2^i;i = i+1;
        try
            solvertime = tic;
            if options.verbose;
                fprintf([' (ok), ' num2str(upper*the_sign)]);
            end
            [sol, flag] = P{upper};
            diagnostic.solvertime = diagnostic.solvertime + toc(solvertime);
        catch
            upper
            diagnostic.problem = -1;
            diagnostic.info = yalmiperror(diagnostic.problem,'BISECTION');
            return
        end
        if upper > 1e6
            if options.verbose;
                     fprintf([' (oj). Giving up, always feasible!\n']);
            end
            diagnostic.problem = 2;
            diagnostic.info = yalmiperror(diagnostic.problem,'BISECTION');
            return
        end
    end
    if options.verbose;
        fprintf([' (fail).']);
    end
end

if options.verbose
    fprintf(['\n']);
end

% Perform bisection
iter = 1;
if options.verbose
    disp('Iteration  Lower bound    Test           Upper bound    Gap          Solver status at test');
end
while upper - lower > options.bisection.absgaptol
    test = (upper + lower)/2;
    solvertime = tic;
    [sol, flag] = P{test};
    diagnostic.solvertime = diagnostic.solvertime + toc(solvertime);
    if options.verbose;
        if options.bisection.switchedsign
            L = -upper;
            T = -test;
            U = -lower;
        else
            L = lower;
            T = test;
            U = upper;
        end
        if flag
           if flag~=1 && ~any(isnan(sol))
                assign(x,sol);assign(Objective,test);res = check(Constraints);
                if min(res) >= -1e-6
                  fprintf(' %4.0f :   %12.5E   %12.5E   %12.5E   %12.5E  %s\n',iter,L,T,U,U-L,[yalmiperror(flag) '(looks ok)'] );
                  flag = 0;
                else
                  fprintf(' %4.0f :   %12.5E   %12.5E   %12.5E   %12.5E  %s\n',iter,L,T,U,U-L,[yalmiperror(flag) '(assumed infeasible)']);
                end
           else
            fprintf(' %4.0f :   %12.5E   %12.5E   %12.5E   %12.5E  %s\n',iter,L,T, U,U-L,yalmiperror(flag));
           end
        else
            fprintf(' %4.0f :   %12.5E   %12.5E   %12.5E   %12.5E  %s\n',iter,L,T, U,U-L,yalmiperror(flag));
        end
    end
    if flag == 0
       working_sol = sol;
       optimal = test;
       lower = test;
    else
        upper = test;
    end
    iter = iter + 1;
end
if options.bisection.switchedsign
    optimal = -optimal;
end
if isempty(working_sol)
    diagnostic.problem = 1;
else
    % Assign computed solution
    assign(x,working_sol);
    assign(Objective,optimal);
end
if options.verbose
    if diagnostic.problem==0
        disp(['Bisection terminated successfully with objective ' num2str(optimal)]);
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function plot_bootfitloglike_a2(mycat,time,timef,bootloops,maepi)
    % Plots Ncum (observed vs. modeled) for specified time windows (choose model for the learning period)
    %
    % plot_bootfitloglike_a2(mycat,time,timef,bootloops,ZG.maepi);
    % --------------------------------------------------
    % Plots Ncum observed vs. Ncum modeled for specified time windows
    % with choosing the model for the learning period
    % Input variables:
    % mycat         : earthquake catalog
    % time      : learning period fo fit Omori parameters [days]
    % timef     : forecast period [days]
    % bootloops : Number of bootstraps
    % ZG.maepi     : mainshock
    %
    % J.Woessner, S. Wiemer

%FIXME mCat ZG.maepi still treated as arrays
report_this_filefun();
    % Surpress warnings from fmincon
    % warning off;

    if ~ensure_mainshock()
        return
    end
    date_matlab = datenum(mycat.Date);
    date_main = datenum(ZG.maepi.Date);
    time_aftershock = date_matlab-date_main;

% Select biggest aftershock earliest in time, but more than 1 day after mainshock
    fDay = 1; %days
    vSel = (mycat.Date > ZG.maepi.Date + days(fDay)) & mycat.Date<= ZG.maepi.Date+days(time);
    mCat = mycat.subset(vSel);
    vSel = mCat.Magnitude == max(mCat.Magnitude);
    vBigAf = mCat(vSel,:);
    if sum(vSel) > 1
        vBigAf.sort('Date')
        vBigAf = vBigAf.subset(1);
    end
    date_biga = datenum(vBigAf.Date);
    fT1 = date_biga - date_main; % Time of big aftershock

    % Aftershock times
    l = time_aftershock(:) > 0;
    tas = time_aftershock(l);
    eqcatalogue = mycat.subset(l);

    % time_as: Learning period
    l = tas <= time;
    time_as=tas(l);

    % Times up to the forecast time
    lf = tas <= time+timef ;
    time_asf= [tas(lf) ];
    time_asf=sort(time_asf);

    % % Calculate p,c,k for dataset
    % prompt  = {'Enter model number (1:pck, 2:pckk, 3:ppckk, 4:ppcckk:'};
    % title   = 'Model selection for fitting aftershock sequence';
    % lines= 1;
    % def     = {'1'};
    % answer  = inputdlg(prompt,title,lines,def);
    % nMod = str2double(answer{1});
    % Calculate fits of different models
    % see OmoriModel for enumaration details

    MOL_models = sort(enumeration('OmoriModel'));

    mRes = [];
    % Modified Omori law (pck)

    for nMod=1:numel(MOL_models) % do this for each model
        [pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL] = bruteforceloglike_a2(time_as, fT1, MOL_models(nMod));
        mRes(nMod,:) = [nMod, pval1, pval2, cval1, cval2, kval1, kval2, fAIC, fL];
    end

    % Select best fitting model by AIC
    vSel = (mRes(:,8)==min(mRes(:,8)));
    mRes = mRes(vSel,:);
    if length(mRes(:,1)) > 1
        vSel = (mRes(:,1)==min(mRes(:,1)));
        mRes = mRes(vSel,:);
    end
    % Model to use for bootstrapping as of lowest AIC to observed data
    nMod = OmoriModel(mRes(1,1));
    pval1= mRes(1,2);
    pval2= mRes(1,3);
    cval1= mRes(1,4);
    cval2= mRes(1,5);
    kval1= mRes(1,6);
    kval2= mRes(1,7);

    % Calculate goodness of fit with KS-Test and RMS
    [H,P,KSSTAT,fRMS] = calc_llkstest_a2(time_as,fT1,pval1, pval2, cval1, cval2, kval1, kval2, nMod);

    % Calculate uncertainty and mean values of p,c,and k
    [mMedModF, mStdL, loopout] = brutebootloglike_a2(time_as, time_asf, bootloops,fT1,nMod);
    pmed1 = mMedModF(1,1);
    pmed2 = mMedModF(1,3);
    cmed1 = mMedModF(1,5);
    cmed2 = mMedModF(1,7);
    kmed1 = mMedModF(1,9);
    kmed2 = mMedModF(1,11);

    % Start plotting
    if (~isnan(pval1) && ~isnan(pval2))

        figure_w_normalized_uicontrolunits('Numbertitle','off','Name','Forecast aftershock occurence')
        loopout(:,end+1) = 0; % add a column

        pval1s = loopout(:,1);
        pval2s = loopout(:,2);
        cval1s = loopout(:,3);
        cval2s = loopout(:,4);
        kval1s = loopout(:,5);
        kval2s = loopout(:,6);

        cumnr_model = OmoriModel.doForecast(nMod, time_asf, pval1s, cval1s, kval1s, fT1, kval2s, pval2s, cval2s);

        loopout(:,9)=max(cumnr_model);

        pfloop = plot(time_asf,cumnr_model,'color',[0.8 0.8 0.8]);
        set(gca,'NextPlot','add')

        warning('This uses only the last values')
        % 2nd moment of bootstrap number of forecasted number of events
        fStdBst = std(loopout(:,9),1,'omitnan');
        %
        % Plot the forecast ...
        cumnrf = (1:length(time_asf))';

        cumnr_modelf = OmoriModel.doForecast(nMod, time_asf, pval1s(end), cval1s(end), kval1s(end), fT1, kval2s(end), pval2s(end), cval2s(end));

        time_asf=sort(time_asf);
        cumnr_modelf=sort(cumnr_modelf);

        pf1 =  plot(time_asf,cumnr_modelf,'g-.','Linewidth',2);
        set(gca,'NextPlot','add')
        %pf2 =  plot(time_asf,cumnrf, 'b-','Linewidth',2);
        %
        % Plot the  fit to the observed data
        % Cumulative number of observed events
        cumnr = (1:length(time_as))';
        cumnr_model = [];
        cumnr_model= OmoriModel.doForecast(nMod, time_as, pval1s(end), cval1s(end), kval1s(end), fT1, kval2s(end), pval2s(end), cval2s(end));

        time_as=sort(time_as);
        cumnr_model=sort(cumnr_model);
        p1 = plot(time_as,cumnr_model,'r','Linewidth',2,'Linestyle','--');
        set(gca,'NextPlot','add');
        p2 = plot(time_as,cumnr,'b','Linewidth',2,'Linestyle','--');

        % Plot the forecast from median value
        cumnr_modelmed = [];
        cumnr_modelmed= OmoriModel.doForecast(nMod, time_asf, pmed1, cmed1, kmed1, fT1, kmed2, pmed2, cmed2);

        time_asf=sort(time_asf);
        cumnr_modelmed=sort(cumnr_modelmed);
        pmedmod =  plot(time_asf,cumnr_modelmed,'y-.','Linewidth',2);

        % Plot observed events in forecast period from endpoint of modeled events in learning period
        vSel = time_asf >= max(time_as);
        vCumnr_forecast = cumnrf(vSel,:);
        vTime_forecast = time_asf(vSel,:);
        % Difference of modelled and observed number of events at time_as
        fDiff_timeas = cumnr_modelmed(length(time_as))-cumnrf(length(time_as));
        vCumnr_forecast = vCumnr_forecast+fDiff_timeas;
        pf3 = plot(vTime_forecast, vCumnr_forecast,'m-.','Linewidth',2);

        xlabel('Time [days]')
        ylabel('Cumulative number of aftershocks')
        xlim([0 max(time_asf)]);

        % Plot standard deviation from bootstrap
        ps2=errorbar(max(time_asf),max(cumnr_modelmed),fStdBst,fStdBst);
        set(ps2,'Linewidth',2,'Color',[1 0 0])
        %
        %     % Title
        % Find amount of events in forecast period for modeled data
        nummod2 = max(cumnr_modelf)-cumnr_modelf(length(time_as));
        nummod = max(cumnr_modelmed)-cumnr_modelmed(length(time_as));
        %     % Find amount of  events in forecast period for observed data
        l = time_asf <=time+timef & time_asf > time;
        numreal = sum(l); % observed number of aftershocks
        %     fRc_Flaw = (numreal-nummod)/sigma;
        fRc_Bst2 = (numreal-nummod2)/fStdBst;
        fRc_Bst = (numreal-nummod)/fStdBst;

        % Round values for output
        pval1 = round(pval1,2);
        pval2 = round(pval2,2);
        cval1 = round(cval1,3);
        cval2 = round(cval2,3);
        kval1 = round(kval1,1);
        kval2 = round(kval2,1);
        pmed1 = round(pmed1,2); mStdL(1,1) = round(mStdL(1,1),2);
        pmed2 = round(pmed2,2); mStdL(1,2) = round(mStdL(1,2),2);
        cmed1 = round(cmed1,3); mStdL(1,3) = round(mStdL(1,3),3);
        cmed2 = round(cmed2,3); mStdL(1,4) = round(mStdL(1,4),3);
        kmed1 = round(kmed1,1); mStdL(1,5) = round(mStdL(1,5),2);
        kmed2 = round(kmed2,1); mStdL(1,6)= round(mStdL(1,6),2);
        fAIC = round(fAIC,2);
        fStdBst = round(fStdBst,2);
        fRc_Bst = round(fRc_Bst,2);

        % Set line for learning period
        yy = get(gca,'ylim');
        plot([max(time_as) max(time_as)],[0 yy(2)],'k-.')
        textX = max(time_asf)*0.1;
        textY = @(z) yy(2)* z;

        switch nMod
            case OmoriModel.pck
                string1 = sprintf('p = %g; c = %g; k = %g', pval1, cval1, kval1);
                string3 = sprintf('pm = %g+-%g; cm = %g+-%g; km = %g+-%g', mped1, mStdL(1,1), cmed1, mStdL(1,3), kmed1, mStdL(1,5));
                % string3=['pm = ' num2str(pmed1) '+-' num2str(mStdL(1,1)) '; cm = ' num2str(cmed1) '+-' num2str(mStdL(1,3)) '; km = ' num2str(kmed1) '+-' num2str(mStdL(1,5))];
                text(textX,textY(0.9),string1,'FontSize',8);
                text(textX,textY(0.8),string3,'FontSize',8);
            case  OmoriModel.pckk
                string1 = sprintf('p = %g; c = %g; k1 = %g; k2 = %g', pval1, cval1, kval1, kval2);
                string3 = sprintf('pm = %g+-%g; cm = %g+-%g; km1 = %g+-%g; km2 = %g+-%g',...
                    mped1, mStdL(1,1), cmed1, mStdL(1,3), kmed1, mStdL(1,5), kmed2, mStdL(1,6));
                text(textX,textY(0.9),string1,'FontSize',8);
                text(textX,textY(0.8),string3,'FontSize',8);
            case OmoriModel.ppckk
                string1 = sprintf('p1 = %g; c = %g; k1 = %g', pval1, cval1, kval1);
                string2 = sprintf('p2 = %g; k2 = %g', pval2, kval2);
                string3 = sprintf('pm1 = %g+-%g; cm = %g+-%g; km1 = %g+-%g', mped1, mStdL(1,1), cmed1, mStdL(1,3), kmed1, mStdL(1,5));
                string4 = sprintf('pm2 = %g+-%g; km2 = %g+-%g', mped2, mStdL(1,2), kmed2, mStdL(1,6));
                text(textX,textY(0.9),string1,'FontSize',8);
                text(textX,textY(0.85),string2,'FontSize',8);
                text(textX,textY(0.8),string3,'FontSize',8);
                text(textX,textY(0.75),string4,'FontSize',8);
            case OmoriModel.ppcckk
                string1 = sprintf('p1 = %g; c1 = %g; k1 = %g', pval1, cval1, kval1);
                string2 = sprintf('p2 = %g; c2 = %g; k2 = %g', pval2, cval2, kval2);
                string3 = sprintf('pm1 = %g+-%g; cm1 = %g+-%g; km1 = %g+-%g', mped1, mStdL(1,1), cmed1, mStdL(1,3), kmed1, mStdL(1,5));
                string4 = sprintf('pm2 = %g+-%g; cm2 = %g+-%g; km2 = %g+-%g', mped2, mStdL(1,2), cmed2, mStdL(1,4), kmed2, mStdL(1,6));
                text(textX,textY(0.9),string1,'FontSize',8);
                text(textX,textY(0.85),string2,'FontSize',8);
                text(textX,textY(0.8),string3,'FontSize',8);
                text(textX,textY(0.75),string4,'FontSize',8);
        end
        string=sprintf('\sigma(Bst) = %g Rc = %g',fStdBst, fRc_Bst);%' Rc(Obfit) = ' num2str(fRc_Bst2)];
        text(textX,yy(2)*0.1,string,'FontSize',8);

        sAIC = sprintf('AIC = %g',fAIC);
        text(textX,yy(2)*0.05,sAIC,'FontSize',8);

        paf = plot(fT1, 0,'h','MarkerFaceColor',[1 1 0],'MarkerSize',12,'MarkerEdgeColor',[0 0 0] );
        sGoodfit = sprintf('KS Test: H = %g, KS statistic = %g P value = %g; RMS = %g', H, KSSTAG, P,fRMS);
        text(textX,yy(2)*0.15,sGoodfit,'FontSize',8);
        % Legend
        legend([p2 p1 pf1 pf3 pmedmod min(ps2) paf],'data','model to data','forecast','observed', 'Mean Bst-model', '\sigma (Bst)','Sec. AF', 'location', 'Best');
    else
        disp('no result')
    end

author: CelsoReyes


text:
function [allUniqNames,allNames] = TCGAfilename2patient(allNames,dashPosition)

% dashPosition is 3 for patient-level and 6 for slide level

% extract patient ID from TCGA format string
for i=1:numel(allNames)
    currName = allNames{i};
    dashes = strfind(currName,'-');
    if numel(dashes)<dashPosition
       % add the last point (empirically tested 29 Jan 2019)
       points = strfind(currName,'.');
       dashes = [dashes,points(end)];
    end
    allNames{i} = currName(1:(dashes(dashPosition)-1));
end

allUniqNames = unique(allNames);

end

author: jnkather

repo: MSIfromHE

save_path: github-repos/MATLAB/jnkather-MSIfromHE


text:
%addpath(genpath('lib'));
%addpath(genpath('IO'));

addpath(genpath('lib/awf'));
addpath(genpath('lib/nyu'));
addpath(genpath('lib/visualindex'));
addpath(genpath('lib/estimateRigidTransform'));
addpath(genpath('lib/peter'));
addpath(genpath('lib/kdtree'));
addpath(genpath('lib/sift'));
addpath(genpath('lib/arcball'));
addpath(genpath('lib/icp'));

warning('off', 'images:imshow:magnificationMustBeFitForDockedFigure');

%fprintf('test\n');
dbstop if error
%clc;

run('lib/vlfeat/toolbox/vl_setup');

author: jianxiongxiao

repo: ProfXkit


text:
%% DEMO_febio_0073_deformable_cylinders_contact_01
% Below is a demonstration for:
%
% * Building geometry for a slab with hexahedral elements, and a
% triangulated sphere.
% * 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
% * indentation
% * contact, sliding, sticky, friction
% * rigid body constraints
% * hexahedral elements, hex8
% * triangular elements, tri3
% * slab, block, rectangular
% * sphere
% * static, solid
% * hyperelastic, Ogden
% * displacement logfile
% * stress logfile

%%

clear; close all; clc;

%% Plot settings
fontSize=15;
faceAlpha1=0.8;
faceAlpha2=0.3;
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 displacement
febioLogFileName_force=[febioFebFileNamePart,'_force_out.txt']; %Log file name for exporting force
febioLogFileName_stress=[febioFebFileNamePart,'_stress_out.txt']; %Log file name for exporting stress

%Geometry parameters
cylRadius1=0.5;
cylRadius2=1;
cylLength=1;
pointSpacing1=(cylRadius1*2*pi)/10*[1 1];
pointSpacing2=(cylRadius1*2*pi)/10*[1 1];

appliedForce=-0.0025;

%Material parameter set
E_youngs1=0.1; %Material Young's modulus
nu1=0.35; %Material Poisson's ratio

E_youngs2=0.1; %Material Young's modulus
nu2=0.35; %Material Poisson's ratio

materialDensity=1e-9;

% FEA control settings
tTotal=1;
numTimeSteps=50; %Number of time steps desired
max_refs=50; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=15; %Optimum number of iterations
max_retries=5; %Maximum number of retires
dtmin=(tTotal/numTimeSteps)/100; %Minimum time step size
dtmax=(tTotal/numTimeSteps); %Maximum time step size
symmetric_stiffness=0;
analysisType='STATIC';%'DYNAMIC';

%Contact parameters
contactPenalty=10;
laugon=0;
minaug=1;
maxaug=10;
fric_coeff=0.01;

%% Creating model geometry and mesh

[meshOutput1]=hexMeshCylinder(cylRadius1,cylLength,pointSpacing1);
E1=meshOutput1.elements;
V1=meshOutput1.nodes;
F1=meshOutput1.faces;
Fb1=meshOutput1.facesBoundary;
Cb1=meshOutput1.boundaryMarker;

[meshOutput2]=hexMeshCylinder(cylRadius2,cylLength,pointSpacing2);
E2=meshOutput2.elements;
V2=meshOutput2.nodes;
F2=meshOutput2.faces;
Fb2=meshOutput2.facesBoundary;
Cb2=meshOutput2.boundaryMarker;

R=euler2DCM([0.5*pi 0 0]);

V1=V1*R;
V1(:,3)=V1(:,3)+cylRadius1;

V2=V2*R;
V2(:,3)=V2(:,3)-cylRadius2;

%% Join node sets

V=[V1;V2]; %Join nodes
E2=E2+size(V1,1); %Shift indices for second set
F2=F2+size(V1,1); %Shift indices for second set
Fb2=Fb2+size(V1,1); %Shift indices for second set
Fb=[Fb1;Fb2];
Cb=[Cb1;Cb2+1+max(Cb1(:))];
Nb1=patchNormal(Fb1,V);
Nb2=patchNormal(Fb2,V);
F=[F1;F2];
E=[E1;E2];

[F,V,~,indFix]=mergeVertices(F,V);
E=indFix(E);
E1=indFix(E1);
E2=indFix(E2);
Fb1=indFix(Fb1);
Fb2=indFix(Fb2);
Fb=indFix(Fb);

%%
% Visualization

cFigure; hold on;

gpatch(Fb1,V,Cb1,'k',faceAlpha1);
patchNormPlot(Fb1,V);

gpatch(Fb2,V,Cb2,'k',faceAlpha1);
patchNormPlot(Fb2,V);

axisGeom(gca,fontSize); camlight headlight;
colormap(gjet(250)); icolorbar;
gdrawnow;

%% Define contact surfaces

nz=[0 0 1];
logicDown1=dot(Nb1,nz(ones(size(Nb1,1),1),:),2)<0;
logicUp2=dot(Nb2,nz(ones(size(Nb2,1),1),:),2)>0;

F_contact_primary=Fb2(Cb2==0 & logicUp2,:) ;
F_contact_secondary=Fb1(Cb1==0 & logicDown1,:);

% Plotting surface models
cFigure; hold on;
title('Contact sets and normal directions','FontSize',fontSize);

gpatch(Fb,V,'kw','none',faceAlpha2);
hl(1)=gpatch(F_contact_secondary,V,'g','k',1);
patchNormPlot(F_contact_secondary,V);
hl(2)=gpatch(F_contact_primary,V,'b','k',1);
patchNormPlot(F_contact_primary,V);

legend(hl,{'Secondary','Primary'});

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

%% Define boundary conditions

VFb1=patchCentre(Fb1,V);
VFb2=patchCentre(Fb2,V);

%Supported nodes
logicPrescribe1=VFb1(:,3)>=cylRadius1 & Cb1==0;
logicSupport2=VFb2(:,3)<=-cylRadius2 & Cb2==0;

bcPrescribeList=unique(Fb1(logicPrescribe1,:));
bcSupportList=unique(Fb2(logicSupport2,:));

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

gpatch(Fb,V,Cb,'none',faceAlpha2);

hl2(1)=plotV(V(bcSupportList,:),'k.','MarkerSize',markerSize);
hl2(2)=plotV(V(bcPrescribeList,:),'r.','MarkerSize',markerSize);

legend(hl2,{'BC support','BC prescribe'});

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=analysisType;
febio_spec.Control.time_steps=numTimeSteps;
febio_spec.Control.step_size=tTotal/numTimeSteps;
febio_spec.Control.solver.max_refs=max_refs;
febio_spec.Control.solver.max_ups=max_ups;
febio_spec.Control.solver.symmetric_stiffness=symmetric_stiffness;
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='neo-Hookean';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.E=E_youngs1;
febio_spec.Material.material{1}.v=nu1;
febio_spec.Material.material{1}.density=materialDensity;

materialName2='Material2';
febio_spec.Material.material{2}.ATTR.name=materialName2;
febio_spec.Material.material{2}.ATTR.type='neo-Hookean';
febio_spec.Material.material{2}.ATTR.id=2;
febio_spec.Material.material{2}.E=E_youngs2;
febio_spec.Material.material{2}.v=nu2;
febio_spec.Material.material{2}.density=materialDensity;

%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(E1,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=E1; %The element matrix

partName2='Part2';
febio_spec.Mesh.Elements{2}.ATTR.name=partName2; %Name of this part
febio_spec.Mesh.Elements{2}.ATTR.type='hex8'; %Element type
febio_spec.Mesh.Elements{2}.elem.ATTR.id=size(E1,1)+(1:1:size(E2,1))'; %Element id's
febio_spec.Mesh.Elements{2}.elem.VAL=E2; %The element matrix

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

nodeSetName2='bcPrescribeList';
febio_spec.Mesh.NodeSet{2}.ATTR.name=nodeSetName2;
febio_spec.Mesh.NodeSet{2}.node.ATTR.id=bcPrescribeList(:);

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

febio_spec.MeshDomains.SolidDomain{2}.ATTR.name=partName2;
febio_spec.MeshDomains.SolidDomain{2}.ATTR.mat=materialName2;

% -> Surfaces
surfaceName1='contactSurface1';
febio_spec.Mesh.Surface{1}.ATTR.name=surfaceName1;
febio_spec.Mesh.Surface{1}.quad4.ATTR.id=(1:1:size(F_contact_secondary,1))';
febio_spec.Mesh.Surface{1}.quad4.VAL=F_contact_secondary;

surfaceName2='contactSurface2';
febio_spec.Mesh.Surface{2}.ATTR.name=surfaceName2;
febio_spec.Mesh.Surface{2}.quad4.ATTR.id=(1:1:size(F_contact_primary,1))';
febio_spec.Mesh.Surface{2}.quad4.VAL=F_contact_primary;

% -> Surface pairs
contactPairName='Contact1';
febio_spec.Mesh.SurfacePair{1}.ATTR.name=contactPairName;
febio_spec.Mesh.SurfacePair{1}.primary=surfaceName2;
febio_spec.Mesh.SurfacePair{1}.secondary=surfaceName1;

%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';

febio_spec.Boundary.bc{2}.ATTR.type='fix';
febio_spec.Boundary.bc{2}.ATTR.node_set=nodeSetName2;
febio_spec.Boundary.bc{2}.dofs='x,y';

% febio_spec.Boundary.bc{3}.ATTR.type='prescribe';
% febio_spec.Boundary.bc{3}.ATTR.node_set=nodeSetName2;
% febio_spec.Boundary.bc{3}.dof='z';
% febio_spec.Boundary.bc{3}.scale.ATTR.lc=1;
% febio_spec.Boundary.bc{3}.scale.VAL=-1;
% febio_spec.Boundary.bc{3}.relative=0;

febio_spec.Loads.nodal_load{1}.ATTR.name='PrescribedForceZ';
febio_spec.Loads.nodal_load{1}.ATTR.type='nodal_load';
febio_spec.Loads.nodal_load{1}.ATTR.node_set=nodeSetName2;
febio_spec.Loads.nodal_load{1}.dof='z';
febio_spec.Loads.nodal_load{1}.scale.ATTR.lc=1;
febio_spec.Loads.nodal_load{1}.scale.VAL=appliedForce/numel(bcPrescribeList);

%Contact section
febio_spec.Contact.contact{1}.ATTR.type='sliding-elastic';
febio_spec.Contact.contact{1}.ATTR.surface_pair=contactPairName;
febio_spec.Contact.contact{1}.two_pass=1;
febio_spec.Contact.contact{1}.laugon=laugon;
febio_spec.Contact.contact{1}.tolerance=0.2;
febio_spec.Contact.contact{1}.gaptol=0;
febio_spec.Contact.contact{1}.minaug=minaug;
febio_spec.Contact.contact{1}.maxaug=maxaug;
febio_spec.Contact.contact{1}.search_tol=0.01;
febio_spec.Contact.contact{1}.search_radius=0.01*sqrt(sum((max(V,[],1)-min(V,[],1)).^2,2));
febio_spec.Contact.contact{1}.symmetric_stiffness=0;
febio_spec.Contact.contact{1}.auto_penalty=1;
febio_spec.Contact.contact{1}.penalty=contactPenalty;
febio_spec.Contact.contact{1}.fric_coeff=fric_coeff;

%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; tTotal 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=',';

febio_spec.Output.logfile.node_data{2}.ATTR.file=febioLogFileName_force;
febio_spec.Output.logfile.node_data{2}.ATTR.data='Rx;Ry;Rz';
febio_spec.Output.logfile.node_data{2}.ATTR.delim=',';

febio_spec.Output.logfile.element_data{1}.ATTR.file=febioLogFileName_stress;
febio_spec.Output.logfile.element_data{1}.ATTR.data='s3';
febio_spec.Output.logfile.element_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='internal';%'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)]);

    %%
    % Importing element stress from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_stress),1,1);

    %Access data
    E_stress_mat=dataStruct.data;

    %%

    % Importing element stress from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_force),1,1);

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

    [CV]=faceToVertexMeasure(E,V,E_stress_mat(:,:,end));

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    title('$\sigma_{3}$ [MPa]','Interpreter','Latex')
    hp=gpatch(Fb,V_DEF(:,:,end),CV,'k',0.8); %Add graphics object to animate
    hp.Marker='.';
    hp.MarkerSize=markerSize2;
    hp.FaceColor='interp';

    axisGeom(gca,fontSize);
    colormap(flipud(turbo(250))); colorbar;
    caxis([min(E_stress_mat(:)) max(E_stress_mat(:))]);
    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

        [CV]=faceToVertexMeasure(E,V,E_stress_mat(:,:,qt));

        %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),CV}; %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:
function [ x, seed ] = uniform_sample ( a, b, seed )

%*****************************************************************************80
%
%% UNIFORM_SAMPLE samples the Uniform PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    A < B.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X, a sample of the PDF.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  [ cdf, seed ] = r8_uniform_01 ( seed );

  x = uniform_cdf_inv ( cdf, a, b );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function Umap=FM2map(im,U,H)
% Unpack fuzzy-membership funcitons to produce membership maps.
%
% INPUT ARGUMENTS:
%   - im    : N-dimensional grayscale image in integer format.
%   - U     : L-by-c array of fuzzy class memberships, where c is the
%             number of classes and L is the intensity range of the input
%             image, such that L=numel(min(im(:)):max(im(:))). See
%             FastFCMeans for more info.
%   - H     : image histogram returned by FastFCMeans function.
%
% OUTPUT:
%   - Umap  : membership maps. Umap has the same size as the input image
%             plus an additional dimension to account for c classes. For
%             example, if im is a 2D M-by-N image then U will be
%             M-by-N-by-c array where U(:,:,i) is a membership map for the
%             i-th class.
%
% AUTHOR    : Anton Semechko (a.semechko@gmail.com)
% DATE      : May.2013
%

if nargin<3 || isempty(H)

    % Intensity range
    Imin=double(min(im(:)));
    Imax=double(max(im(:)));
    I=(Imin:Imax)';

    % Intensity histogram
    H=hist(double(im(:)),I);
    H=H(:);

end

% Unpack memberships
Umap=zeros(sum(H),size(U,2));
i1=1; i2=0;
for i=1:numel(H)
    i2=i2+H(i);
    Umap(i1:i2,:)=repmat(U(i,:),[H(i) 1]);
    i1=i2+1;
end

% Find the positional mapping
[~,idx]=sort(im(:), 'ascend');
[~,idx]=sort(idx(:),'ascend');

% Reshape membership maps to match image dimensions.
Umap=reshape(Umap(idx,:),[size(im) size(U,2)]);

author: Sable