text:
function lung_image = PTKFillCoronalHoles(lung_image, is_right, reporting)
    % PTKFillCoronalHoles. Operates on each coronal slice, applying a closing
    % filter then filling interior holes
    %
    %
    %
    %     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.
    %

    reporting.ShowProgress('Filling holes');

    max_non_coronal_voxel_size = max(lung_image.VoxelSize(2:3));

    opening_size = 10/max_non_coronal_voxel_size;
    closing_size = 10/max_non_coronal_voxel_size;

    lung_image.AddBorder(10);
    for coronal_index = 11 : lung_image.ImageSize(1) - 10
        reporting.UpdateProgressStage(coronal_index - 11, lung_image.ImageSize(1) - 20);
        coronal_slice = lung_image.GetSlice(coronal_index, PTKImageOrientation.Coronal);
        if ~isempty(is_right)
            coronal_slice = OpenOrClose(coronal_slice, is_right, opening_size, closing_size, reporting);
        end
        coronal_slice = imclose(coronal_slice, strel('disk', round(closing_size)));
        coronal_slice = imfill(coronal_slice, 'holes');
        lung_image.ReplaceImageSlice(coronal_slice, coronal_index, PTKImageOrientation.Coronal);
    end
    lung_image.RemoveBorder(10);
    reporting.CompleteProgress;
end

function mask = OpenOrClose(mask, is_right, opening_size, closing_size, reporting)
    closed_image = imclose(mask, strel('disk', round(closing_size)));
    opened_image = imopen(mask, strel('disk', round(opening_size)));

    threshold = PTKImage(~(closed_image & opened_image));
    threshold.AddBorder(1);
    image_size = threshold.ImageSize;

    raw_image = zeros(image_size);
    if (is_right)
        raw_image(2, 2:end-1, 2) = 2;
        raw_image(end-1, 2:end-1, 2) = 1;
    else
        raw_image(2, 2:end-1, 2) = 1;
        raw_image(end-1, 2:end-1, 2) = 2;
    end
    raw_image(2:end-1, end-1, 2) = 2;
    raw_image(2:end-1, 2, 2) = 2;
    right_border_indices = threshold.LocalToGlobalIndices(find(raw_image == 1));
    other_border_indices = threshold.LocalToGlobalIndices(find(raw_image == 2));

    start_points = {right_border_indices, other_border_indices};

    reporting.PushProgress;
    regions = PTKMultipleRegionGrowing(threshold, start_points, reporting);
    reporting.PopProgress;

    regions.RemoveBorder(1);
    closed_region = (regions.RawImage == 1);
    mask(closed_region) = closed_image(closed_region);
    mask(~closed_region) = opened_image(~closed_region);
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [vITRS,rotMat]=ITRS2TIRS(vTIRS,TT1,TT2,xpyp)
%%ITRS2TIRS Rotate vectors from the International terrestrial Reference
%            System (ITRS) into the  Terrestrial Intermediate Reference
%            System (TIRS). The ITRS is essentially the WGS-84 coordinate
%            system: it defines locations with respect to the crust of a
%            non-rotating Earth, where the z axis passes through a fixed
%            point on the surface. On the other hand, the TIRS is nearly
%            the same except the z axis is the axis of rotation of the
%            Earth, which slowly varies over time. Note that the velocity
%            conversion does not include the (small) centrifugal effect of
%            polar motion.
%
%INPUTS: x The NXnumVec collection of vectors in TIRS coordinates to
%          convert (units do not matter). N can be 3, or 6. If the vectors
%          are 3D, then they are position. 6D vectors are assumed to be
%          position and velocity. Since the TIRS and ITRS co-rotate, there
%          is no Coriolis effect to add. Also, the accelerations due to the
%          wobble of the rotation axis over time are not considered. These
%          accelerations are very small. Thus, the function just rotates
%          both halves of the vector.
% TT1, TT2 Two parts of a Julian date given in terrestrial time (TT).
%          The units of the date are days. The full date is the sum of
%          both terms. The date is broken into two parts to provide more
%          bits of precision. It does not matter how the date is split.
%     xpyp xpyp=[xp;yp] are the polar motion coordinates in radians
%          including the effects of tides and librations. If this
%          parameter is omitted or if an empty matrix is passed, the value
%          from the function getEOP will be used.
%
%OUTPUTS: vITRS The NXnumVec vector of values of x rotated from the ITRS
%               into the TIRS.
%        rotMat The 3X3 rotation matrix used to rotate vectors from the
%               ITRS into the TIRS.
%
%The conversion functions from the International Astronomical Union's
%(IAU) Standard's of Fundamental Astronomy library are put together to get
%the necessary rotation matrix.
%
%The algorithm can be compiled for use in Matlab  using the
%CompileCLibraries function.
%
%The algorithm is run in Matlab using the command format
%[vITRS,rotMat]=ITRS2TIRS(vTIRS,TT1,TT2)
%or if more parameters are known, using the format
%[vITRS,rotMat]=ITRS2TIRS(vTIRS,TT1,TT2,xpyp);
%
%Different celestial coordinate systems are compared in [1].
%
%REFERENCES:
%[1] D. F. Crouse, "An Overview of Major Terrestrial, Celestial, and
%    Temporal Coordinate Systems for Target Tracking," Formal Report,
%    Naval Research Laboratory, no. NRL/FR/5344--16-10,279, 10 Aug. 2016,
%    173 pages.
%
%March 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

error('This function is only implemented as a mexed C or C++ function. Please run CompileCLibraries.m to compile the function for use.')

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [out] = recharge_6(p1,p2,S,dt)
%recharge_6

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Recharge to fulfil evaporation demand if the receiving
%               store is below a threshold
% Constraints:  f <= S/dt
%               S >= 0      prevents complex numbers
% @(Inputs):    p1   - time coefficient [d-1]
%               p2   - non-linear scaling [mm]
%               S    - current storage [mm]
%               dt   - time step size [d]

out = min(max(S/dt,0),p1.*max(S,0).^p2);

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function value = dasum ( n, x, incx )

%*****************************************************************************80
%
%% DASUM takes the sum of the absolute values of a vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 June 2005
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for Fortran Usage,
%    Algorithm 539,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, September 1979, pages 308-323.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, real X(*), the vector to be examined.
%
%    Input, integer INCX, the increment between successive entries of X.
%    INCX must not be negative.
%
%    Output, real VALUE, the sum of the absolute values of X.
%
  value = sum ( abs ( x(1:incx:1+(n-1)*incx) ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [h, axh]=barplot_colored(data,varargin)
% Make a barplot of data with error bars, with colors specified by colormap
% or color string.
%
% :Usage:
% ::
%
%    [h, axh]=barplot_colored(data, [optional arguments])
%
% this is a good function though
% within-subject error bars now added; use 'within'
% ::
%
%    [bar_handles, axis_handle]=barplot_colored(data,varargin)
%
% :Input arguments: Optional
%
%   **'within':**
%        Do within-subject STE bars, average obs x variable interaction
%        Loftus and Masson 1994 style.
%
%   Strings, followed by values for each:
%
%   :COLOR CONTROL:
%
%   **'colormap':**
%        followed by colormap name to use
%
%   **'colors':**
%        followed by cell array of colors per bar; supercedes colormap
%
%   :Display items:
%     - 'fontsize'
%     - 'title'
%     - 'XTickLabels'
%     - 'ylabel'
%     - 'xlabel'
%
%   :Bar Locations:
%     - 'x' : followed by x values for bars (locations)
%
% NOTE: For this function, keywords must be even-numbered argument entries,
% e.g., arg 2, 4, 6.  Odd argument entries are values.
%
% For example: This works, and you need the extra empty arg after 'within'
% ::
%
%    [h1, s1] = barplot_colored(pexp1, 'within', ' ', 'title', 'Pattern expression', 'XTickLabels', dat.Y_names, 'x', 1:nterms);
%
% You can assign arbitrary colors to bars by setting the colormap:
% ::
%
%    [h, s] = barplot_colored([corr_temp corr_rep]);
%    cm = [1 .5 0; .5 0 1];
%    colormap(cm)
%
% Example: A grouped barplot
% ::
%
%    dat = rand(20, 4);
%    create_figure('bars');
%    [h1, s1] = barplot_colored(dat, 'x', [1 2 4 5]);
%    % set(h2, 'BarWidth', .9)
%
% Change colormap:
% ::
%
%    [h1, s1] = barplot_colored(dat, 'x', [1 2 4 5], 'colormap', 'summer');
%
% Enter values:
% ::
%
%    colors = {[.8 .25 .25] [.8 .5 .25] [.4 .5 .8] [.25 .25 .9]};
%    [h1, s1] = barplot_colored(dat, 'x', [1 2 4 5], 'colors', colors);
%
% Set X Tick Label:
% ::
%
%    [h1, s1] = barplot_colored(dat, 'XTicklabels', {'A' 'B' 'C' 'D'});
%
% See also: barplot_columns, lineplot_columns

if iscell(data)
    for k=1:length(data)
        means(k)=nanmean(data{k});
        stderr(k)=nanstd(data{k})/sqrt(sum(~isnan(data{k})));
    end
else
    means=mean(data);
    stderr=std(data)/sqrt(size(data,1));
end

x = 1:length(means); % can replace x values
myfontsize = 18;
mytitle = '';

for k=1:length(varargin)
    if strcmp(varargin{k},'colormap')

        eval(['colormapfun=@' varargin{k+1} ';'])
    end
    if ischar(varargin{k}) && strcmp(varargin{k},'title')
        mytitle=varargin{k+1};
    end
    if ischar(varargin{k}) && strcmp(lower(varargin{k}),'fontsize')
        myfontsize=varargin{k+1};
    end
    if ischar(varargin{k}) && (strcmp(lower(varargin{k}),'xticklabel') || strcmp(lower(varargin{k}),'xticklabels'))
        XTickLabel=varargin{k+1};
    end
    if ischar(varargin{k}) && strcmp(lower(varargin{k}),'ylabel')
        Ylabel=varargin{k+1};
    end
    if ischar(varargin{k}) && strcmp(lower(varargin{k}),'xlabel')
        Xlabel=varargin{k+1};

    end
    if ischar(varargin{k}) && strcmp(varargin{k},'x')
        x=varargin{k+1};
    end
    if ischar(varargin{k}) && strcmp(varargin{k},'colors')

        colors=varargin{k+1};

    end

    if ischar(varargin{k}) && strcmp(varargin{k},'within')
        if iscell(data), error('Within error bars not implemented for cell input data'); end

                stderr = barplot_get_within_ste(data);
                stderr = repmat(stderr, 1, length(means));

    end
end

k = size(x, 2);

if ~exist('colormapfun','var')
    colormapfun=@hsv;
end

if ~exist('colors','var') || isempty(colors)
    colors = colormapfun(k);
end

if ~iscell(colors)
    for i = 1:size(colors, 1)
        tmp{i} = colors(i, :);
    end
    colors = tmp;
end

axh = gca;
set(axh, 'FontSize', myfontsize)
title(mytitle);

for i = 1:k
    hold on
    h(i) = bar(x(i), means(i));
    set(h(i), 'FaceColor', colors{i});
end

% This is old code, Matlab has changed...2014/5 update
% h=bar(x, means);
% s=get(h,'Children');
%
% colormap(colormapfun(length(means)));
% set(s,'CData',1:length(means));

errorbar(x, means,stderr,'k','LineWidth',2,'LineStyle','none')
set(gca,'Xlim',[0 max(x)+1])
if exist('XTickLabel','var')
    set(gca,'XTickLabel',XTickLabel, 'XTick', x)
else
    set(gca,'XTickLabel',[])
end
if exist('mytitle','var')
    title(mytitle)
end
if exist('Ylabel','var')
    ylabel(Ylabel)
end
if exist('Xlabel','var')
    xlabel(Xlabel)
end

hold off

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function Save_EWT2D_LP(ewt2d,name)

%=========================================================================
%
% function Save_EWT2D_LP(ewt2d,name)
%
% Save each subband image in separate files named 'LPxname.png' where n
% is the subband number.
% The images are renormalized with respect to the min and max of all
% images.
%
% Author: Jerome Gilles
% Institution: UCLA - Department of Mathematics
% Year: 2013
% Version: 1.0
% ========================================================================

minI=255.0;
maxI=0.0;

% Find the min and max over all images
for n=1:length(ewt2d)
    minI=min(minI,min(ewt2d{n}(:)));
    maxI=max(maxI,max(ewt2d{n}(:)));
end

% Plot each subband with the same normalization
for n=1:length(ewt2d)
%        gname=sprintf('LP%d%s.png',n,name);
        gname=sprintf('R%d%s.png',n,name);

        imwrite((ewt2d{n}-min(ewt2d{n}(:)))/(max(ewt2d{n}(:))-min(ewt2d{n}(:))),gname,'png');
%        imwrite((ewt2d{n}-minI)/(maxI-minI),gname,'png');
end

author: Sable

repo: mcbench-benchmarks


text:
function styleLabelAsWarn(hLabel,warnStr)
%styleLabelAsWarn Summary of this function goes here
%   Detailed explanation goes here
    lblStr = strWrapForLabel(hLabel, warnStr);

    set(hLabel,'Visible','on');
    set(hLabel,'BackgroundColor',[1,1,0]);
    set(hLabel,'ForegroundColor',[0,0,0]);
    set(hLabel,'String',lblStr);
    set(hLabel,'TooltipString','');
    set(hLabel,'FontWeight','bold');

    setappdata(hLabel,'tooltipstring',warnStr);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function results = eq(A, B)
   % when using wildcards, use ismember
   switch class(B)
      case 'scnlobject'
         results = [A.tag] == [B.tag];
      case 'ChannelTag'
         results = [A.tag] == [B];
      otherwise
         error('scnlobject:eq:UnrecognizedComparison',...
            'Don''t know how to compare scnlobject with %s\n',class(B));
   end
end

author: geoscience-community-codes

repo: GISMO


text:
% Coupled nonlinear PDE's
% Buckley Leverett equation
% dependent variables: pressure and water saturation
% Prepared for educational purposes by ** AAE **
% TVD scheme, really sharp fronts! Not too bad computationally.
% Written by Ali A. Eftekhari
% Last checked: June 2021
clc
%% define the geometry
Nx = 200; % number of cells in x direction
Ny = 150; % number of cells in y direction
W = 50; % [m] length of the domain in x direction
H = 30; % [m] length of the domain in y direction
m = createMesh2D(Nx, Ny, W, H); % creates a 2D mesh
%% define the physical parametrs
p0 = 1e5; % [bar] pressure
pin = 50e5; % [bar] injection pressure at the left boundary
sw0 = 0; % initial water saturation
swin = 1;
mu_oil = 10e-3; % [Pa.s] oil viscosity
mu_water = 1e-3; % [Pa.s] water viscosity
% reservoir
k0 = 2e-12; % [m^2] average reservoir permeability
phi0 = 0.2; % average porosity
clx=0.05;
cly=0.05;
V_dp=0.6; % Dykstra-Parsons coef.
perm_val= field2d(Nx,Ny,k0,V_dp,clx,cly);
k=createCellVariable(m, perm_val);
phi=createCellVariable(m, phi0);
krw0 = 1;
kro0 = 1;
nw = 2;
no = 2;
krw = @(sw)(krw0*sw.^nw);
dkrwdsw = @(sw)(krw0*nw*sw.^(nw-1));
kro = @(sw)(kro0*(1-sw).^no);
dkrodsw = @(sw)(-kro0*no*(1-sw).^(no-1));
lw = geometricMean(k)/mu_water;
lo = geometricMean(k)/mu_oil;
Lw = @(sw)(krw(sw));
Lo = @(sw)(k/mu_oil*kro(sw));
dLwdsw = @(sw)(k/mu_water*dkrwdsw(sw));
dLodsw = @(sw)(k/mu_oil*dkrodsw(sw));
%% Define the boundaries
BCp = createBC(m); % Neumann BC for pressure
BCs = createBC(m); % Neumann BC for saturation
% change the left and right boandary to constant pressure (Dirichlet)
BCp.left.a(:)=0; BCp.left.b(:)=1; BCp.left.c(:)=pin;
BCp.right.a(:)=0; BCp.right.b(:)=1; BCp.right.c(:)=p0;
% change the left boundary to constant saturation (Dirichlet)
BCs.left.a(:)=0; BCs.left.b(:)=1; BCs.left.c(:)=1;
%% define the time step and solver properties
dt = 1000; % [s] time step
t_end = 1000*dt; % [s] final time
eps_p = 1e-5; % pressure accuracy
eps_sw = 1e-5; % saturation accuracy
%% define the variables
sw_old = createCellVariable(m, sw0, BCs);
p_old = createCellVariable(m, p0, BCp);
sw = sw_old;
p = p_old;
uw = -gradientTerm(p_old); % an estimation of the water velocity
%% start the main loop
% generate intial pressure profile (necessary to initialize the fully
% implicit solver)
FL = fluxLimiter('SUPERBEE');
t = 0;
while (t<t_end)
    error_p = 1e5;
    error_sw = 1e5;
    % estimation loop: sequential
%     for i = 1:2
%         pgrad = gradientTerm(m, p.value);
%         sw_face = upwindMean(m, -pgrad, sw.value); % average value of water saturation
%         labdao = lo.*funceval(kro, sw_face);
%         labdaw = lw.*funceval(krw, sw_face);
%         labda = labdao+labdaw;
%         Mdiffp1 = diffusionTerm(m, -labda);
%         [Mbcp, RHSbcp] = boundaryCondition(m, BCp);
%         p.value = reshapeCell(m, full((Mdiffp1+Mbcp)\RHSbcp));
%         for j = 1:2
%             pgrad = gradientTerm(m, p.value);
%             sw_face = upwindMean(m, -pgrad, sw.value);
%             labdaw = lw.*funceval(krw, sw_face);
%             dlabdawdsw = lw.*funceval(dkrwdsw, sw_face);
%             [Mbcsw, RHSbcsw] = boundaryCondition(m, BCs);
%             [Mtranssw2, RHStrans2] = transientTerm(m, phi, dt, sw);
%             Mconvsw2 = convectionUpwindTerm(m, -dlabdawdsw.*pgrad);
%             RHSs = RHSbcsw+RHStrans2+divergenceTerm(m, (labdaw-dlabdawdsw.*sw_face).*pgrad);
%             Ms = Mbcsw+Mconvsw2+Mtranssw2;
%             sw.value = reshapeCell(m, full(Ms\RHSs));
%         end
%     end
    % Implicit loop
    while ((error_p>eps_p) || (error_sw>eps_sw))
        % calculate parameters
        pgrad = gradientTerm(p);

        % sw_face = upwindMean(m, -pgrad, sw.value); % average value of water saturation
        sw_face = tvdMean(sw, -pgrad, FL); % average value of water saturation
        labdao = lo.*funceval(kro, sw_face);
        labdaw = lw.*funceval(krw, sw_face);
        dlabdaodsw = lo.*funceval(dkrodsw, sw_face);
        dlabdawdsw = lw.*funceval(dkrwdsw, sw_face);
        labda = labdao+labdaw;
        dlabdadsw = dlabdaodsw+dlabdawdsw;
        % compute [Jacobian] matrices
        [Mconvsw1, RHSconvsw1] = convectionTvdTerm(-dlabdadsw.*pgrad, sw, FL);
%         sw_face = tvdMean(m, sw.value, -pgrad, FL); % average value of water saturation
%         labdao = lo.*funceval(kro, sw_face);
%         labdaw = lw.*funceval(krw, sw_face);
%         dlabdaodsw = lo.*funceval(dkrodsw, sw_face);
%         dlabdawdsw = lw.*funceval(dkrwdsw, sw_face);
%         labda = labdao+labdaw;
%         dlabdadsw = dlabdaodsw+dlabdawdsw;
        [Mconvsw2, RHSconvsw2] = convectionTvdTerm(-dlabdawdsw.*pgrad, sw, FL);

        [Mtranssw2, RHStrans2] = transientTerm(sw_old, dt, phi);
        Mdiffp1 = diffusionTerm(-labda);
        Mdiffp2 = diffusionTerm(-labdaw);
        % Compute RHS values
        RHS1 = divergenceTerm(-dlabdadsw.*sw_face.*pgrad);
        RHS2 = divergenceTerm(-dlabdawdsw.*sw_face.*pgrad);
        % include boundary conditions
        [Mbcp, RHSbcp] = boundaryCondition(BCp);
        [Mbcsw, RHSbcsw] = boundaryCondition(BCs);
        % Couple the equations; BC goes into the block on the main diagonal
        M = [Mdiffp1+Mbcp Mconvsw1; Mdiffp2 Mconvsw2+Mtranssw2+Mbcsw];
        RHS = [RHS1+RHSbcp+RHSconvsw1; RHS2+RHStrans2+RHSbcsw+RHSconvsw2];
        % solve the linear system of equations
        x = M\RHS;
        % separate the variables from the solution
        p_new = reshapeCell(m,full(x(1:(Nx+2)*(Ny+2))));
        sw_new = reshapeCell(m,full(x((Nx+2)*(Ny+2)+1:end)));
        % calculate error values
        error_p = max(max(abs((p_new-p.value)./p_new)));
        error_sw = max(max(abs(sw_new-sw.value)));
        % assign new values of p and sw
        p.value = p_new;
        sw.value = sw_new;
    end
    t=t+dt;
    p_old = p;
    sw_old = sw;
    figure(1);visualizeCells(sw); drawnow;%shading interp
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function kern = matern32KernParamInit(kern)

% MATERN32KERNPARAMINIT MATERN32 kernel parameter initialisation.
% The Matern class of kernels is a wide class with different
% degrees of freedom parameters. This is the specific case where nu
% = 3/2.
%
% Given
% r = sqrt((x_i - x_j)'*(x_i - x_j))
%
% We have
% k(x_i, x_j) = sigma2*(1+sqrt(3)*r/l)*exp(-sqrt(3)*r/l)
%
% The parameters are sigma2, the process variance (kern.variance),
% and l, the length scale (kern.lengthScale).
% FORMAT
% DESC initialises the matern kernel with nu=3/2
%  kernel structure with some default parameters.
% ARG kern : the kernel structure which requires initialisation.
% RETURN kern : the kernel structure with the default parameters placed in.
%
% SEEALSO : kernCreate, kernParamInit
%
% COPYRIGHT : Neil D. Lawrence, 2006

% KERN

kern.variance = 1;
kern.lengthScale = 1;
kern.nParams = 2;

kern.transforms.index = [1 2];
kern.transforms.type = optimiDefaultConstraint('positive');

kern.isStationary = true;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function outn=eq(x,y)

precAndSize

for ii=1:max(ex,ey)
 imag=false;
 if ex==1
  [xrval,xival]=getVals(x,1);
  [yrval,yival]=getVals(y,ii);
 elseif ey==1
  [xrval,xival]=getVals(x,ii);
  [yrval,yival]=getVals(y,1);
 else
  [xrval,xival]=getVals(x,ii);
  [yrval,yival]=getVals(y,ii);
 end
  outn(ii)=mpfr_eq(precision,xrval,yrval);
  if outn(ii)==0, outn(ii)=1; else, outn(ii)=0; end
end % for ii=1:max(ex,

author: opencobra

repo: cobratoolbox


text:
% computes bootstrap bias corrected accelerated bootstrap confidence
% interval and returns if it includes 0 as a logical vector sig
% mu_est - estimated parameter from complete sample
% bs_est - bootstrap estimates
% jk_est - jacknife estimates
% alpha - desired FPR
%
% Written by Bogdan Petre, 2020
function [sig, CI] = BCa(mu_est, bs_est, jk_est, alpha)
    fprintf('Estimating BCa CIs...\n')
    bias = sum(bs_est < mu_est)./length(bs_est);
    z0 = icdf('norm', bias, 0, 1);

    jk_mean = mean(jk_est);
    num = sum((jk_mean - jk_est).^3);
    den = sum((jk_mean - jk_est).^2);
    a = num ./ (6*den.^(3/2));

    zL = z0 + icdf('norm',alpha/2,0,1);
    alpha1 = normcdf(z0 + zL./(1-a.*zL));
    zU = z0 + icdf('norm',1-alpha/2,0,1);
    alpha2 = normcdf(z0 + zU./(1-a.*zL));

    CI = quantile(bs_est, [alpha1, alpha2]);

    sig = prod(sign(CI),2) > 0;
end


text:
function compObjectHypot( aid )
%COMPOBJECTHYPOT Generate object hypothesis
%   You should have run compRoomHypot first

global config;
bufname = config.bufname;

rot_file_name = [bufname config.IMGLIST(aid).name config.rotateImageSuffix];
rotImg = imread(rot_file_name);

% selective search
SS = getSelectiveSearch(rotImg);
AllCuboid = regionBasedHypothesisFromSS(SS, 1024, 2048);
rectangle = rectDetectionFlexFromSS(SS, 1024, 2048);
AllCuboid_NEW = pruneCuboid(rotImg, AllCuboid);
rectangle_NEW = pruneRectangle(rotImg, rectangle);
% rectangle hypothesis
[ rectangle, rectCuboid2, rectCuboid3 ] = rectangleBasedHypothesis( rotImg );
% segmentation hypothesis
[ regionCuboid ] = regionBasedHypothesis( rotImg );
% merge up everything
regionCuboid.views = [regionCuboid.views; AllCuboid_NEW.views];
regionCuboid.xyzBox = [regionCuboid.xyzBox; AllCuboid_NEW.xyzBox];
regionCuboid.score = [regionCuboid.score; AllCuboid_NEW.score];
regionCuboid.count = regionCuboid.count + AllCuboid_NEW.count;

for i = 1:6
    rectangle(i).highPixelBox = [rectangle(i).highPixelBox; rectangle_NEW(i).highPixelBox];
    rectangle(i).xyzBox = [rectangle(i).xyzBox; rectangle_NEW(i).xyzBox];
    rectangle(i).score = [rectangle(i).score; rectangle_NEW(i).score];
end

if config.UPDATE_TO_DISK
    parsave([bufname config.objHypoFileUpdate], ...
        'rectangle', rectangle, 'rectCuboid2', rectCuboid2, 'rectCuboid3', rectCuboid3, ...
        'regionCuboid', regionCuboid);
end

[ rect_score ] = getSegScrOfRect( rotImg, rectangle );
CAN_POOL = generateObjHyps( aid, 20000, rect_score);
if config.UPDATE_TO_DISK
    save([bufname config.rectangleScore], 'rect_score');
    save([bufname config.objHypoFile3D],'CAN_POOL');
end
fprintf('finish.\n');
end

author: zouchuhang

repo: LayoutNet


text:
function MV = ffmInteractionsSparseHF(X,Xbox,Y,Ybox,V,Ibox,green,k,edg,tol)
%+========================================================================+
%|                                                                        |
%|         OPENFFM - LIBRARY FOR FAST AND FREE MEMORY CONVOLUTION         |
%|           openFfm is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2019.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : ffmInteractionsSparseHF.m                     |
%|    #    |   VERSION    : 0.6                                           |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 05.09.2019                                    |
%| ( === ) |   SYNOPSIS   : Sparse product for high-frequency compressible|
%|  `---'  |                leaves                                        |
%+========================================================================+

% Initialisation du produit Matrice-Vecteur
MV = zeros(size(X,1),1,class(V));

% Quadrature spherique de Geggenbauer
[Xq,Wq,l] = ffmQuadratureHF(k,edg,tol);

% Unicite des vecteurs de translation
XY        = Ybox.ctr(Ibox(:,2),:) - Xbox.ctr(Ibox(:,1),:);
[~,Il,It] = unique(floor(XY*1e6),'rows','stable');
Nt        = length(Il);

% Fonctions de transfert
Txy = cell(Nt,1);
for i = 1:Nt
    Txy{i} = ffmTransfertHF(Xq,Wq,XY(Il(i),:),k,l);
end

% Convolution en Y
Vy = cell(size(Ybox.ind));
for i = unique(Ibox(:,2)')
    % Boite centree en Y
    iy = Ybox.ind{i};
    ny = length(iy);
    Ym = Y(iy,:) - ones(ny,1)*Ybox.ctr(i,:);

    % Transformee de Fourier inverse (Ym->Xq)
    Vy{i} = ffmInterpHF(Xq,Ym,V(iy),-1,tol);

    % Derivation du noyau en Y
    if strcmp(green(1:end-1),'grady[exp(ikr)/r]')
        j     = str2double(green(end));
        Vy{i} = -1i*Xq(:,j) .* Vy{i};
    end
end

% Translations
TVy = cell(size(Xbox.ind));
for i = 1:length(TVy)
    TVy{i} = 0;
end
for i = 1:size(Ibox,1)
    TVy{Ibox(i,1)} = TVy{Ibox(i,1)} + Txy{It(i)}.*Vy{Ibox(i,2)};
end

% Convolutions en X
for i = unique(Ibox(:,1)')
    % Boite centree en X
    ix = Xbox.ind{i};
    nx = length(ix);
    Xm = X(ix,:) - ones(nx,1)*Xbox.ctr(i,:);

    % Derivation du noyau en X
    if strcmp(green(1:end-1),'gradx[exp(ikr)/r]')
        j      = str2double(green(end));
        TVy{i} = 1i*Xq(:,j) .* TVy{i};
    end

    % Transformee de Fourier (Xq->X)
    MV(ix) = ffmInterpHF(Xm,Xq,TVy{i},+1,tol);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Xq,Wq,l] = ffmQuadratureHF(k,edg,tol)
% Ordre harmonique (E. Darve)
l = floor( abs(k)*sqrt(3)*edg - log(tol) );

% Quadrature de Gauss sur [-1,1]
u     = 1:l;
u     = u ./ sqrt(4*u.^2 - 1);
[V,x] = eig( diag(u,-1) + diag(u,+1) );
[x,I] = sort(diag(x));
w     = 2*V(1,I)'.^2;

% Quadrature de Gauss sur [0,1] par transformation lineaire
a = 0;   b = 1;
x = 0.5*(b-a)*x + 0.5*(a+b);
w = 0.5*(b-a)*w;

% Quadrature de Gauss en elevation (phi)
phi = acos(2*x(end:-1:1)-1) - pi/2;
Wp  = 2*w(end:-1:1)';

% Quadrature reguliere en azimut (theta)
Ntheta = 2*(l+1);
theta  = 2*pi/Ntheta*(0:Ntheta-1)';
Wt     = 2*pi/Ntheta*ones(Ntheta,1);

% Quadrature spherique par produit tensoriel
[theta,phi] = ndgrid(theta,phi);
theta       = theta(:);
phi         = phi(:);
Wq          = Wt * Wp;
Wq          = Wq(:);

% Coordonnees carthesienne et produit par le nombre d'onde
[x,y,z] = sph2cart(theta,phi,1);
Xq      = k*[x,y,z];

% Securite
if (l>3)
    % Harmoniques spheriques jusqu'a l'ordre 3
    n   = 1;
    Ylm = zeros(length(theta),16);
    for ll = 0:3
        Plm = sqrt((2*ll+1)/(4*pi)) * legendre(ll,cos(pi/2-phi),'sch').';
        for m = -ll:ll
            if m<0
                Ylm(:,n) = Plm(:,abs(m)+1) .* sin(abs(m).*theta);
            else
                Ylm(:,n) = Plm(:,abs(m)+1) .* cos(abs(m).*theta);
            end
            n = n+1;
        end
    end

    % Test de l'integration spherique des harmoniques (produit scalaire)
    [m,n] = size(Ylm);
    if norm( (Ylm' * spdiags(Wq,0,m,m) * Ylm) - eye(n) , 'inf') > 1e-5
        error('ffmQuadratureHF.m - error 1');
    end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function T = ffmTransfertHF(Xq,ws,xy,k,l)
% Fonctions de hankel spherique 1ere espece
besselhs = @(n,z) sqrt(pi./(2.*z)) .* besselh(n+0.5,1,z);

% Distance
kr = k*norm(xy);

% cosinus(angle incidence)
x = Xq*(-xy/kr)';

% Initialisation de la recursion sur (besselhs(kr) * Pl(cos(phi))
P0 = ones(size(x));
P1 = x;
T  = besselhs(0,kr).*P0 + 3i*besselhs(1,kr).*P1;

% Construction recursive des polynomes de Legendre
for ll = 2:l
    P2 = 1/ll .* ( (2*(ll-1)+1).*x.*P1 - (ll-1).*P0 );
    T  = T + (2*ll+1) * 1i^ll * besselhs(ll,kr) .* P2;
    P0 = P1;
    P1 = P2;
end

% Operateur de Translation
T = 1i*abs(k)/(4*pi) .* ws .* T;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function MV = ffmInterpHF(X,Y,V,iflag,tol)
% Dimensions
Nx   = size(X,1);
Ny   = size(Y,1);
Nmax = max(Nx,Ny);
Nmin = min(Nx,Ny);

% DFT non uniforme
if (Nmin<100) || (Nmin<log(Nmax))
    MV = exp(1i*iflag*X*Y') * V;

% FFT non uniforme
else
    % Conversion de type
    if isa(X,'single') || isa(Y,'single') || isa(V,'single')
        sgl = 1;
        X = double(X); Y = double(Y); V = double(V);
    else
        sgl = 0;
    end

    % Initialisation produit Matrice-Vecteur
    MV = zeros(Nx,1) + 1i*zeros(Nx,1);

    % Prevention de coplanarite en X
    Nx         = Nx + 1;
    X(end+1,:) = X(end,:) + tol;
    MV(end+1)  = 0;

    % Prevention de coplanarite en X
    Ny         = Ny + 1;
    Y(end+1,:) = Y(end,:) + tol;
    V(end+1)   = 0;

    % Transformee de Fourier
    ier    = 0;
    mex_id = 'nufft3d3f90(i int[x], i double[], i double[], i double[], i dcomplex[], i int[x], i double[x], i int[x], i double[], i double[], i double[], io dcomplex[], io int[x])';
    MV     = nufft3d(mex_id,Ny,Y(:,1),Y(:,2),Y(:,3),V, iflag, tol, ...
        Nx,X(:,1),X(:,2),X(:,3),MV, ier, 1, 1, 1, 1, 1);

    % Supression dernier point
    MV = MV(1:end-1);

    % Conversion de type
    if sgl
        MV = single(MV);
    end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function isgn = i4vec_compare ( n, a1, a2 )

%*****************************************************************************80
%
%% I4VEC_COMPARE compares two I4VEC's.
%
%  Discussion:
%
%    The lexicographic ordering is used.
%
%  Example:
%
%    Input:
%
%      A1 = ( 2, 6, 2 )
%      A2 = ( 2, 8, 12 )
%
%    Output:
%
%      ISGN = -1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vectors.
%
%    Input, integer A1(N), A2(N), the vectors to be compared.
%
%    Output, integer ISGN, the results of the comparison:
%    -1, A1 < A2,
%     0, A1 = A2,
%    +1, A2 < A1.
%
  isgn = 0;

  k = 1;

  while ( k <= n )

    if ( a1(k) < a2(k) )
      isgn = -1;
      return
    elseif ( a2(k) < a1(k) )
      isgn = + 1;
      return
    end

    k = k + 1;

  end

  return
end


text:
function plane = fit_depth_plane(x_hom, depth)
%FIT_DEPTH_PLANE  Least squares fitting of a plane in the 3D scene projected on
%the image from a set of input points.
%
%   INPUTS:
%
%   -|x_hom|: |P|-by-3 matrix holding homogeneous image-plane coordinates of
%   input points.
%
%   -|depth|: |P|-by-1 vector holding depth values of input points.
%
%   OUTPUTS:
%
%   -|plane|: 3-by-1 vector holding the parameters of the fitted plane, so that
%   for each point [x, y, d] it holds d = [x, y, 1] * plane.

x_hom_T = x_hom.';
plane = (x_hom_T * x_hom) \ (x_hom_T * depth);

end

author: sakaridis


text:
function A = sllinreg(X, Y, varargin)
%SLLINREG Performs Multivariate Linear Regression and Ridge Regression
%
% $ Syntax $
%   - A = sllinreg(X, Y, ...)
%
% $ Arguments $
%   - X:        The matrix of x samples
%   - Y:        The matrix of y samples
%   - A:        The solved transform matrix
%
% $ Description $
%   - A = sllinreg(X, Y, varargin) solves the linear regression problem to
%     get the linear transforms A: (y = Ax + e). The solution is given by
%     the following optimization problem:
%       A = argmin_{A} sum_i ||y_i - A x_i||^2 + sum_k lambda_k ||a_k||^2
%     Here, Y is a dy x n matrix with the i-th column giving y_i;
%           X is a dx x n matrix with the i-th column giving x_i
%           A is a dy x dx transform matrix
%           a_k is the k-row vector of A, which corresponds to the k-th
%                  component of y
%           lambda_k is the regularization weight in ridge regression
%     According to mathematical analysis, the solution is given as
%       A = (Y * X^T) * (X * X^T + diag(lambda))^{-1}
%     You can specify the following properties to control the regression:
%     \*
%     \t    Table  Linear Regression Properties
%     \h        name       &             description
%              'lambdas'   &  The regularization coefficients in ridge
%                             linear regression. If all components share
%                             the same value, then lambdas can be a scalar.
%                             If different values are for different
%                             dimensions, then lambdas can be a dy x 1
%                             column vector. (default = 0)
%              'weights'   &  The sample weights, can be either [] or
%                             an 1 x n row vector.
%              'invparams' &  The parameters for invoking slinvcov to
%                             compute inverse matrix, in the form of
%                             {method, ...}. default = [], means directly
%                             using inv to do inverse.
%     \*
%
% $ Remarks $
%   - To solve the linear problem like y = Ax + b, you have two ways:
%     (1) centralize x and y respectively, and invoke sllinreg, then set
%         b = my - A * mx
%     (2) invoke sllinrega, which is specially designed for such an
%         augmented problem.
%
% $ History $
%   - Created by Dahua Lin, on Sep 15, 2006
%

%% parse and verify input arguments

if nargin < 2
    raise_lackinput('sllinreg', 2);
end

if ~isnumeric(X) || ~isnumeric(Y) || ndims(X) ~= 2 || ndims(Y) ~= 2
    error('sltoolbox:invalidarg', ...
        'The X and Y should be both 2D numeric matrices');
end

[dx, n] = size(X);
[dy, ny] = size(Y);
if n ~= ny
    error('sltoolbox:sizmismatch', ...
        'X and Y contain different numbers of samples');
end

opts.lambdas = 0;
opts.rv = 1e-3;
opts.weights = [];
opts.invparams = [];
opts = slparseprops(opts, varargin{:});

lambdas = opts.lambdas;
if ~isscalar(lambdas) && ~isequal(size(lambdas), [dy, 1])
    error('lambdas should be either a scalar or a dy x 1 column vector');
end

w = opts.weights;
if ~isempty(w)
    if ~isequal(size(w), [1, n])
        error('The sample weights should be a 1 x n row vector');
    end
end

%% main

if isempty(w)
    Xt = X';
else
    Xt = slmulvec(X, w, 2)';
end

M2 = X * Xt;
if ~isequal(lambdas, 0)
    inds = (1:dx)' * (dx+1) - dx;
    M2(inds) = M2(inds) + lambdas;
end

if isempty(opts.invparams)
    invM2 = inv(M2);
else
    invM2 = slinvcov(M2, opts.invparams{:});
end
clear M2;

M1 = Y * Xt;

A = M1 * invM2;

author: lmthang


text:
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function gdpowert
%GDPOWERT Unit test for the function GDPOWER.

%	O. Lemoine - August 1996.

N=256;

% K=1 : constant group delay
[sig,gd]=gdpower(N,1);
t0=N; t=(1:N);
if any(gd~=t0),
 error('gdpower test 1 failed');
end
if any(sig(t~=t0)>sqrt(eps)),
 error('gdpower test 2 failed');
end

% K=2 : linear group delay
[sig,gd]=gdpower(N,2);
ifl=instfreq(sig);
[s,iflaw]=fmlin(N);
if any(abs(ifl(3:N-2)-iflaw(3:N-2))>5e-2),
 error('gdpower test 3 failed');
end

% K=1/2
[sig,gpd]=gdpower(N,1/2);
fnorm=linspace(0.01,0.45,108);
fnorm=fnorm(2:108);
gd=sgrpdlay(sig,fnorm);
if any(abs(gd-gpd(4:110)')/N>5e-2),
 error('gdpower test 4 failed');
end;

% K=0
[sig,gpd]=gdpower(N,0);
fnorm=linspace(0.01,0.45,108);
fnorm=fnorm(2:108);
gd=sgrpdlay(sig,fnorm);
if any(abs(gd-gpd(4:110)')/N>5e-2),
 error('gdpower test 5 failed');
end;

% K=-1
[sig,gpd]=gdpower(N,-1);
fnorm=linspace(0.01,0.45,108);
fnorm=fnorm(2:108);
gd=sgrpdlay(sig,fnorm);
if any(abs(gd-gpd(4:110)')/N>5e-2),
 error('gdpower test 6 failed');
end;

N=221;

% K=1 : constant group delay
[sig,gd]=gdpower(N,1);
t0=N; t=(1:N);
if any(gd~=t0),
 error('gdpower test 7 failed');
end
if any(sig(t~=t0)>sqrt(eps)),
 error('gdpower test 8 failed');
end

% K=2 : linear group delay
[sig,gd]=gdpower(N,2);
ifl=instfreq(sig);
[s,iflaw]=fmlin(N);
if any(abs(ifl(5:N-2)-iflaw(5:N-2))>5e-2),
 error('gdpower test 9 failed');
end

% K=1/2
[sig,gpd]=gdpower(N,1/2);
fnorm=linspace(0.01,0.45,108);
fnorm=fnorm(2:108);
gd=sgrpdlay(sig,fnorm);
if any(abs(gd-gpd(4:110)')/N>5e-2),
 error('gdpower test 10 failed');
end;

% K=0
[sig,gpd]=gdpower(N,0);
fnorm=linspace(0.01,0.45,108);
fnorm=fnorm(2:108);
gd=sgrpdlay(sig,fnorm);
if any(abs(gd-gpd(4:110)')/N>5e-2),
 error('gdpower test 11 failed');
end;

% K=-1
[sig,gpd]=gdpower(N,-1);
fnorm=linspace(0.01,0.45,108);
fnorm=fnorm(2:108);
gd=sgrpdlay(sig,fnorm);
if any(abs(gd-gpd(4:110)')/N>5e-2),
 error('gdpower test 12 failed');
end;

author: HeLiangHIT


text:
function [V, converged, i] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt)
%NEWTONPF  Solves power flow using full Newton's method (power/polar)
%   [V, CONVERGED, I] = NEWTONPF(YBUS, SBUS, V0, REF, PV, PQ, MPOPT)
%
%   Solves for bus voltages using a full Newton-Raphson method, using nodal
%   power balance equations and polar coordinate representation of
%   voltages, given the following inputs:
%       YBUS  - full system admittance matrix (for all buses)
%       SBUS  - handle to function that returns the complex bus power
%               injection vector (for all buses), given the bus voltage
%               magnitude vector (for all buses)
%       V0    - initial vector of complex bus voltages
%       REF   - bus index of reference bus (voltage ang reference & gen slack)
%       PV    - vector of bus indices for PV buses
%       PQ    - vector of bus indices for PQ buses
%       MPOPT - (optional) MATPOWER option struct, used to set the
%               termination tolerance, maximum number of iterations, and
%               output options (see MPOPTION for details).
%
%   The bus voltage vector contains the set point for generator
%   (including ref bus) buses, and the reference angle of the swing
%   bus, as well as an initial guess for remaining magnitudes and
%   angles.
%
%   Returns the final complex voltages, a flag which indicates whether it
%   converged or not, and the number of iterations performed.
%
%   See also RUNPF, NEWTONPF_S_CART, NEWTONPF_I_POLAR, NEWTONPF_I_CART.

%   MATPOWER
%   Copyright (c) 1996-2019, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   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.

%% default arguments
if nargin < 7
    mpopt = mpoption;
end

%% options
tol         = mpopt.pf.tol;
max_it      = mpopt.pf.nr.max_it;
lin_solver  = mpopt.pf.nr.lin_solver;

%% initialize
converged = 0;
i = 0;
V = V0;
Va = angle(V);
Vm = abs(V);

%% set up indexing for updating V
npv = length(pv);
npq = length(pq);
j1 = 1;         j2 = npv;           %% j1:j2 - V angle of pv buses
j3 = j2 + 1;    j4 = j2 + npq;      %% j3:j4 - V angle of pq buses
j5 = j4 + 1;    j6 = j4 + npq;      %% j5:j6 - V mag of pq buses

%% evaluate F(x0)
mis = V .* conj(Ybus * V) - Sbus(Vm);
F = [   real(mis([pv; pq]));
        imag(mis(pq))   ];

%% check tolerance
normF = norm(F, inf);
if mpopt.verbose > 1
    fprintf('\n it    max P & Q mismatch (p.u.)');
    fprintf('\n----  ---------------------------');
    fprintf('\n%3d        %10.3e', i, normF);
end
if normF < tol
    converged = 1;
    if mpopt.verbose > 1
        fprintf('\nConverged!\n');
    end
end

%% attempt to pick fastest linear solver, if not specified
if isempty(lin_solver)
    nx = length(F);
    if nx <= 10 || have_feature('octave')
        lin_solver = '\';       %% default \ operator
    else    %% MATLAB and nx > 10 or Octave and nx > 2000
        lin_solver = 'LU3';     %% LU decomp with 3 output args, AMD ordering
    end
end

%% do Newton iterations
while (~converged && i < max_it)
    %% update iteration counter
    i = i + 1;

    %% evaluate Jacobian
    [dSbus_dVa, dSbus_dVm] = dSbus_dV(Ybus, V);
    [dummy, neg_dSd_dVm] = Sbus(Vm);
    dSbus_dVm = dSbus_dVm - neg_dSd_dVm;

    j11 = real(dSbus_dVa([pv; pq], [pv; pq]));
    j12 = real(dSbus_dVm([pv; pq], pq));
    j21 = imag(dSbus_dVa(pq, [pv; pq]));
    j22 = imag(dSbus_dVm(pq, pq));

    J = [   j11 j12;
            j21 j22;    ];

    %% compute update step
    dx = mplinsolve(J, -F, lin_solver);

    %% update voltage
    if npv
        Va(pv) = Va(pv) + dx(j1:j2);
    end
    if npq
        Va(pq) = Va(pq) + dx(j3:j4);
        Vm(pq) = Vm(pq) + dx(j5:j6);
    end
    V = Vm .* exp(1j * Va);
    Vm = abs(V);            %% update Vm and Va again in case
    Va = angle(V);          %% we wrapped around with a negative Vm

    %% evalute F(x)
    mis = V .* conj(Ybus * V) - Sbus(Vm);
    F = [   real(mis([pv; pq]));
            imag(mis(pq))   ];

    %% check for convergence
    normF = norm(F, inf);
    if mpopt.verbose > 1
        fprintf('\n%3d        %10.3e', i, normF);
    end
    if normF < tol
        converged = 1;
        if mpopt.verbose
            fprintf('\nNewton''s method power flow (power balance, polar) converged in %d iterations.\n', i);
        end
    end
end

if mpopt.verbose
    if ~converged
        fprintf('\nNewton''s method power flow (power balance, polar) did not converge in %d iterations.\n', i);
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [J, grad] = costFunctionReg(theta, X, y, lambda)
%COSTFUNCTIONREG Compute cost and gradient for logistic regression with regularization
%   J = COSTFUNCTIONREG(theta, X, y, lambda) computes the cost of using
%   theta as the parameter for regularized logistic regression and the
%   gradient of the cost w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta

predictions =  sigmoid(X*theta);

leftPart = -y' * log(predictions);

rightPart = (1 - y') * log(1 - predictions);

thetaZero = theta;

thetaZero(1) = 0;

lambaCostPart = (lambda / (2 * m)) * sum(thetaZero .^ 2);

lambdaGradPart = lambda / m * thetaZero;

J = (1 / m) * (leftPart - rightPart) + lambaCostPart;

grad = ((1/m) * (X' * (predictions - y))) + lambdaGradPart;

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

end

author: vugsus

repo: coursera-machine-learning


text:
% *************************************************************************
% Code for setting the numerical values of nonlinear terms
% *************************************************************************
function p = updateonenonlinearbound(p,changed_var)
if ~isempty(p.bilinears)
    impactedVariables = find((p.bilinears(:,2) == changed_var) | (p.bilinears(:,3) == changed_var));
    x = p.bilinears(impactedVariables,2);
    y = p.bilinears(impactedVariables,3);
    z = p.bilinears(impactedVariables,1);
    x_lb = p.lb(x);
    x_ub = p.ub(x);
    y_lb = p.lb(y);
    y_ub = p.ub(y);
    bounds = [x_lb.*y_lb x_lb.*y_ub x_ub.*y_lb x_ub.*y_ub];
    p.lb(z) = max([p.lb(z) min(bounds,[],2)],[],2);
    p.ub(z) = min([p.ub(z) max(bounds,[],2)],[],2)';
    p.lb(impactedVariables(x==y)<0) = 0;
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [sMin,sMax,s0]=projEllipseOntoLine(z,A,gammaVal,x0,v)
%%PROJELLIPSEONTOLINE Determine the orthogonal projection of an ellipse or
%                     ellipsoid onto a line.
%
%INPUTS: z The numDimX1 center of the ellipsoid.
%        A A numDimXnumDim symmetric, positive definite matrix that
%          specifies the size and shape of the ellipse or ellipsoid, where
%          a point zp is on the ellipse/ellipsoid if
%          (zp-z)'*A*(zp-z)=gammaVal.
% gammaVal The threshold for declaring a point to be in the ellipsoid. If
%          an empty matrix is passed, the default value of 1 is used.
%     x0,v Two numDimX1 vectors that specify the line. A point y is on the
%          line if y=x0+s*v, where s is a scalar quantity.
%
%OUTPUTS: sMin,sMax The range of values of s on the line y=x0+s*v onto
%                   which the ellipsoid is orthogonally projected. These
%                   are scalar quantities.
%
%The mathematics for projecting an ellipsoid onto a line are described in
%Section 12 of [1].
%
%EXAMPLE:
% A=[1,0;
%    0,10];
% M=[0.413074198133900,  0.910697373904216;
%    -0.910697373904216,   0.413074198133900];%A rotation matrix
% A=M*A*M';
% z=[5;6];
% x0=[0;0];
% v=[1;2];
% [sMin,sMax,s0]=projEllipseOntoLine(z,A,[],x0,v);
% figure(1)
% clf
% hold on
% axis([-0.5, 7.5, 2, 10])%Same size in x and y
% axis square
% drawEllipse(z,A,1,'b','linewidth',2)
% %Draw a segment of the line.
% xStart=x0+v;
% xEnd=x0+5*v;
% plot([xStart(1);xEnd(1)],[xStart(2);xEnd(2)],'-k','linewidth',4)
% %Draw the segment of the orthogonal projection onto the line.
% xStart=x0+sMin*v;
% xEnd=x0+sMax*v;
% plot([xStart(1);xEnd(1)],[xStart(2);xEnd(2)],'-r','linewidth',2)
% %Draw the line from the center of the ellipse to the place on the line onto
% %which it is orthogonally projected.
% xS0=x0+s0*v;
% plot([z(1);xS0(1)],[z(2);xS0(2)],'--c')
%
%REFERENCES:
%[1] S. B. Pope, "Algorithms for ellipsoids," Cornell University, Tech.
%    Rep. FDA-08-01, Feb. 2008. [Online].
%    Available: https://tcg.mae.cornell.edu/pubs/Pope_FDA_08.pdf
%
%May 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(isempty(gammaVal))
   gammaVal=1;
end

A=A/gammaVal;

L=chol(A,'lower');

s0=v'*(z-x0)/(v'*v);
w=norm(L\v/(v'*v));

sMin=s0+w;
sMax=s0-w;
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory


text:
filename='Chair_Tetrahedra';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'holes';
cost = {'compliance'};
weights = [1];
constraint = {'volume'};
optimizer = 'HAMILTON-JACOBI'; incrementFactor = 1;
filterType = 'P1';
line_search_initiator = 'STANDARD';

HJiter0 = 1;
e2 = 100;
N_holes = [24 11 11];
R_holes = 0.4;
phase_holes = [0 0 0];
nsteps = 5;

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

Vfrac_initial = 0.4;
optimality_initial = 5e-2;
constr_initial = 5e-2;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;


text:
function e = rmModelSearchFit(p,Y,trends,params,rawrss);
% rmModelSearchFit - actual fit function of rmSearchFit
%
% error = rmModelSearchFit(p,Y,trends,params);
%
% Basic fit with several hard limits on where the values can
% go. When encountering a set limit error will go to infinity. We
% may want to consider a more smooth limit.
%
% 2006/06 SOD: wrote it.

% sigma should be > 0 and < sigmaRatioInfVal
% This is a hard border beyond which the estimates cannot go.
if p(3)<=0 | p(3)>params.analysis.sigmaRatioInfVal, e = realmax; return; end;

% don't estimate too far away from the stimulus size, say 10x
if sqrt(p(1).^2+p(2).^2) > 10*max([params.stim(:).stimSize]), e = realmax; return; end;

% input check for 1 or 2 Gaussians
if numel(p) == 3,
  gaussianId = [p(3) p(3) 0 p(1) p(2)];
else,
  % also check that second Gaussian is at least twice as big as the
  % first one
  if p(4)< 2.*p(3) | p(4)>params.analysis.sigmaRatioInfVal, e = realmax; return; end;
  gaussianId = [p(3) p(3) 0 p(1) p(2);p(4) p(4) 0 p(1) p(2)];
end;

% make prediction
[pred, weight] = rfMakePrediction(params,gaussianId);

% Another hard border that depends on the relative overlap with
% stimulus window and penalizes pRFs too far away or too
% large. We'll do this only for the 1st pRF.
if weight(1)<0.01, e = realmax; return; end;

% fit
X = [pred trends];
b = pinv(X)*Y;

% Last sanity checks:
% The first pRF should be positive.
if b(1)<0, e = realmax; return; end;
% For two Gaussians the center should always be positive.
if numel(p) == 4,
  if b(1)+b(2)<=0, e = realmax; return; end;
end;

% RSS weighted by amount within of RF in stimulus window
% e = sum((Y - X*b).^2) ./ weight;
e = norm(Y - X*b).^2;

% if rawrss is given scale rss relative to rawrss,
% otherwise report actuall rss.
if nargin > 4,
  e = e./rawrss.*100;
end;

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [mvector,scan] = bes(RF, pw, scantype, var1, var2, var3, var4, phase, M0, T1, T2)

% Bloch Equation Simulator
%
% Required Arguments
%  RF        4xN definition of radio frequency pulse (phase, amplitude,
%               duration, gradient)
%  pw        duration of pulse in milliseconds
%  scantype  scan type: 'frequency', 'B1max', 'temporal'
%  B1max scan
%     var1     frequency / position (if gradient supplied in RF)
%     var2     upper bound B1
%     var3     lower bound B1
%     var4     steps
%  Frequency scan
%     var1     B1max
%     var2     upper bound frequency / position (if gradient supplied)
%     var3     lower bound frequency / position (if gradient supplied)
%     var4     steps
%  Temporal scan
%     var1     B1max
%     var2     frequency / position (if gradient supplied)
%     var3     ignored
%     var4     ignored
% Optional Arguments
%  phase    phase at which the RF pulse is applied (default 0)
%  M0       three component vector specifying the initial magnetization
%              a 3xSTEPS matrix is used to define a different MO for each
%              steps. Default is [0; 0; 1]
%  T1       T1 decay time (default INF)
%  T2       T2 decay time (default INF)
%
% Output
%  mvector  3xSTEPS matrix specifying the magnetization at each point
%  scan     vector of frequency, B1max, or time steps modelled
%
% Units: B1 (kHz), frequencies (kHz), Phase (degrees), Time (ms),
%        Gradients (G/cm), Position (cm), gyro = 4.25763888 kHz/G

% Author: L. Martyn Klassen
% Copyright 2003 Robarts Research Institute
% This program is copyrighted worldwide by the Robarts Research
% Institute.  Any distribution, copying or redistribution is expressly
% forbidden.
% 03/04/24 Added support for gradient changes during RF

if nargin < 5
   error('BES requires five input arguments.');
end

% Initialize constants
scantype = lower(scantype(1));
grad_flag = 0;

if nargin < 11
   T2 = Inf;
   if nargin < 10
      T1 = Inf;
      if nargin < 9
         M0 = [0; 0; 1];
         if nargin < 8
            phase = 0;
            if (nargin < 7) && (scantype ~= 't')
               error('BES requires seven input arguments for specified scan type.');
            end
         end
      end
   end
end

% If only one column is defined in the RF pulse it is an amplitude
% modulated pulse, and the phase, duration columns have to be added
if ndims(RF) > 2
   error('BES: RF pulse cannot have more than 2 dimensions.');
end

sizRF = size(RF);
if all(sizRF > 4)
   error('BES: RF pulse has too many columns (phase, amplitude, duration, grad)')
elseif sizRF(2) > 4
   RF = RF';
end

switch size(RF, 2)
   case 1
      % Single column is RF magnitude so you have to add phase, duration,
      % and gradients
      RF(:,2) = RF;
      RF(:,1) = phase;
      RF(:,3) = 1;
      RF(:,4) = 0;
   case 2
      % Two columns are phase and magnitude, add duration.
      % Add the desired phase to the waveform:
      RF(:,1) = RF(:,1) + phase;
      RF(:,3) = 1;
   case 3
      % Three columns are phase, magnitude and duration.
      % Add the desired phase to the waveform:
      RF(:,1) = RF(:,1) + phase;
   case 4
      % Convert from G/cm to kHz/cm
      RF(:,4) = RF(:,4) .* 4.25763888;
      grad_flag = 1;
      % Add the desired phase to the waveform:
      RF(:,1) = RF(:,1) + phase;
end

% Convert phase to radians
RF(:,1) = RF(:,1) .* (pi / 180);

% Find out how many cyles are required to step through the RF pulse
cycles = size(RF, 1);

% Scale the RF amplitude to be between -1 and 1
RF(:,2) = RF(:,2)./max(abs(RF(:,2)));

% Kill any points with zero duration
RF = RF(RF(:,3) ~= 0, :);

% The smallest duration allowed is 1.0
minduration = min(RF(:,3));
maxduration = max(RF(:,3));
if (minduration < 1)
   error('BES: Pulse step duration should be zero or positive integer.');
end

% Adjust durations to be integers
RF(:,3) = round(RF(:,3));

% Pulse width has to be a scalar
pw = pw(1);

% Define the smallest time interval used in the pulse
dt = pw./sum(RF(:,3));

% A time scan uses the pulse squence to define the progression
% through time and not steps
if scantype == 't'
   steps = 1;
else
   steps = var4;
end

% Define M0
if isempty(M0)
   % Default is alignment along z axis
   M0 = [0; 0; 1];
   if scantype ~= 't'
      mvector = M0(:,ones(steps, 1));
   end
elseif prod(size(M0)) == 3
   % One vector must be repeated to size of the scan
   % Check to make sure that M0 is not too large
   if sum(M0.^2) > 1.0
      error('BES: Initial magnetization is greater than 1.0')
   else
      M0 = M0(:);
      if scantype ~= 't'
         mvector = M0(:,ones(steps, 1));
      end
   end
elseif size(M0, 2) ~= steps
   if scantype ~= 't'
      % Otherwise M0 must be defined at each point in the scan
      error('BES: Initial M0 vector does not agree with number of steps (3xsteps).')
   else
      % Time requires it to 3x1 or 1x3 - previous if statement
      error('BES: Initial MO for time scan must be 3x1 vector')
   end
else
   if scantype == 't'
      error('BES: Initial MO for time scan must be 3x1 vector')
   end
   if any(sum(M0.^2, 1) > 1.000001)
      error('BES: Initial magnetization is greater than 1.0')
   end
   mvector = M0;
end

% Make sure lower limit is less than upper limit
% switch if they are the wrong order
if var2 < var3
   temp = var2;
   var2 = var3;
   var3 = temp;
   clear temp;
end

if (T2 > T1)
   error('BES: First Law Violation, T2 must be less than T1.');
end

% Calculate the matrices for the T1, T2 decay
if T1 == 0
   if T2 == 0
      A = [0 0 0; 0 0 0; 0 0 1];
   else
      A = [exp(-dt/T2) 0 0; 0 exp(-dt/T2) 0; 0 0 1];
   end
   B = [0; 0; 0];
else
   if T2 == 0
      A = [0 0 0; 0 0 0; 0 0 1-exp(-dt/T1)];
   else
      A = [exp(-dt/T2) 0 0; 0 exp(-dt/T2) 0; 0 0 exp(-dt/T1)];
   end
   B = [0; 0; 1-exp(-dt/T1)];
end

% Setup for the three different types of simulations
switch scantype
case 'f'   % Scan through frequency

   % Fixed value is B1max which has to be scaled for the RF pulse description
   % Also convert it to rad/millisecond
   B1max = var1.*(2*pi);

   % Limits are in frequency
   interval = (var2 - var3)/(steps-1);
   scan = [var3:interval:var2];

   % Rotate the matrix to align 1st RF with X
   cosa = cos(-RF(1,1));
   sina = sin(-RF(1,1));
   mvector = [cosa sina 0; -sina cosa 0; 0 0 1]*mvector;

   % Small tip angle vector
   theta = RF(:,2)*B1max.*dt;
   cosa = cos(theta);
   sina = sin(theta);

   % Free precession angle (rad)
   freq = scan .* (2*pi*dt);
   phi = freq;
   cos2a = cos(phi);
   sin2a = sin(phi);

   % Alignment angles
   dRF = -diff(RF(:,1));
   dRF(end+1) = RF(end, 1);

   % Loop through the RF pulse
   for t = 1:cycles
      % Compute the rotation matrix
      Rx=[1 0 0; 0 cosa(t) sina(t); 0 -sina(t) cosa(t)];

      % Added in the decay matrix
      Rx = A * Rx;

      if grad_flag
         phi = freq .* RF(t,4);
         cos2a = cos(phi);
         sin2a = sin(phi);
      end

      % For durations longer than one, repeat the small angle
      % approximation the prerquisite number of times
      for m = 1:RF(t,3)
         % Using small tip angle approximation to calculate the
         % new magnetization vector

         % Rotate the matrix by theta about X and let decay
         mvector = Rx*mvector;

         % Add the constant decay term
         mvector(3,:) = mvector(3,:) + B(3);

         % Free Precession angle with RF alignment on last increment
         if (m == RF(t,3)) & (dRF(t) ~= 0)
            phi2 = phi + dRF(t);
            cos3a = cos(phi2);
            sin3a = sin(phi2);

            temp = cos3a.*mvector(1,:) + sin3a.*mvector(2,:);
            mvector(2,:) = -sin3a.*mvector(1,:) + cos3a.*mvector(2,:);
            mvector(1,:) = temp;
         else
            % Rotate the vector by phi about Z
            temp = cos2a.*mvector(1,:) + sin2a.*mvector(2,:);
            mvector(2,:) = -sin2a.*mvector(1,:) + cos2a.*mvector(2,:);
            mvector(1,:) = temp;
         end
      end
   end

case 'b'    % Scan through B1max

   % Fixed value is Frequency (rad/ms)
   freq = var1*2*pi;

   % Limits are in B1max which has to be scaled for the RF pulse description
   % and converted to rad/ms.
   interval = (var2 - var3)/(steps-1);
   B1max = [var3:interval:var2].*(2*pi);

   % Rotate the matrix to align 1st RF with X
   cosa = cos(-RF(1,1));
   sina = sin(-RF(1,1));
   mvector = [cosa sina 0; -sina cosa 0; 0 0 1]*mvector;

   % Free Precession rotation angle (rad)
   freq = freq*dt;
   phi = freq;
   % Compute the free precession matrix
   cosa = cos(phi);
   sina = sin(phi);
   Rfree =[cosa sina 0; -sina cosa 0; 0 0 1];

   % Add the decay term as well
   Rfree = A * Rfree;

   % Calculate RF alignment angles
   dRF = -diff(RF(:,1));
   % Last one must realign to zero
   dRF(end+1) = RF(end,1);
   cos2a = cos(phi + dRF);
   sin2a = sin(phi + dRF);

   % Loop through the RF pulse
   for t = 1:cycles
      % Compute theta for given RF power
      theta = B1max*RF(t,2)*dt;
      cosa = cos(theta);
      sina = sin(theta);

      if grad_flag
         phi = freq * RF(t,4);
      end

      % For durations longer than one, repeat the small angle
      % approximation the prerequisite number of times
      for m = 1:RF(t,3)
         % Using small tip angle approximation to calculate the
         % new magnetization vector

         % Rotate the matrix by theta about X
         temp = cosa.*mvector(2,:) + sina.*mvector(3,:);
         mvector(3,:) = -sina.*mvector(2,:) + cosa.*mvector(3,:);
         mvector(2,:) = temp;

         % On last time through apply free precession and rotate the matrix
         % to undo the previous alignment with the RF and also align the next
         % RF with X, otherwise just apply free precession
         if grad_flag
            % Compute the free precession matrix and apply
            cosa = cos(phi + dRF(t));
            sina = sin(phi + dRF(t));
            mvector = (A * [cosa sina 0; -sina cosa 0; 0 0 1])*mvector;
         elseif (m == RF(t,3)) & (dRF(t) ~= 0)
            mvector=(A * [cos2a(t) sin2a(t) 0; -sin2a(t) cos2a(t) 0; 0 0 1])*mvector;
         else
            mvector = Rfree*mvector;
         end
         mvector(3,:) = mvector(3,:) + B(3);
      end
   end

   % Convert scan from rad/ms to kHz
   scan = B1max./(2*pi);

case 't'
   % Scan through time

   % Fixed value are B1max and freq
   B1max = var1.*(2*pi);
   freq = var2.*2*pi;

   mvector = zeros(3, sum(RF(:,3))+1);
   count = 1;

   % Store the starting magnetization
   mvector(:,count) = M0;

   % Free Precession rotation
   freq = freq*dt;
   phi = freq;

   % Loop through the RF pulse
   for t = 1:cycles
      if grad_flag
         phi = freq .* RF(t,4);
      end

      % Rotate the matrix to align RF with X
      cosa = cos(-RF(t,1));
      sina = sin(-RF(t,1));
      Rz1 = [cosa sina 0; -sina cosa 0; 0 0 1];

      % Compute the RF rotation matrix
      theta = B1max*RF(t,2)*dt;
      cosa = cos(theta);
      sina = sin(theta);
      Rx=[1 0 0; 0 cosa sina; 0 -sina cosa];

      % Also rotate the matrix to undo the previous alignment with the RF
      % and free precession
      cosa = cos(phi + RF(t,1));
      sina = sin(phi + RF(t,1));
      Rz2=[cosa sina 0; -sina cosa 0; 0 0 1];

      % Add in the decay terms
      Rz2 = A * Rz2;

      % For durations longer than one, repeat the small angle
      % approximation the prerquisite number of times
      for m = 1:RF(t,3)
         % Using small tip angle approximation to calculate the
         % new magnetization vector

         % RF alignment
         M0 = Rz1*M0;

         % Rotate the matrix by theta about X
         M0 = Rx*M0;

         % Free precesion, realignment and decay
         M0 = Rz2*M0 + B;

         % Store the magnetization vector
         count = count + 1;
         mvector(:,count) = M0;
      end
   end

   scan = [0:sum(RF(:,3))].*dt;

otherwise
   error('BES: Scan type is frequency, B1max, or time');
end

return

author: CIC-methods

repo: FID-A

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


text:
function f = kappa_inf(tau, P)
%  This function computes approximation to kappa in the limit of n->inf
%
x0 = norminv( (tau+1)/2, 0, 1);
VP = 2 * (1 + 2*P);
VQ = 2 * (1 + P)^2;
f = 2*normcdf(sqrt(VP/VQ) * x0, 0, 1) - 1;

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function quadrule_test28 ( )

%*****************************************************************************80
%
%% TEST28 tests LEGENDRE_SET_X0_01 and SUM_SUB.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2009
%
%  Author:
%
%    John Burkardt
%
  order_max = 8;

  nfunc = func_set ( 'COUNT', 'DUMMY' );

  a = 0.0;
  b = 1.0;

  nsub = 1;

  xlo =  0.0;
  xhi =  1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST28\n' );
  fprintf ( 1, '  LEGENDRE_SET_X0_01 sets up Gauss-Legendre quadrature\n' );
  fprintf ( 1, '    for integrating F(X) over [0,1];\n' );
  fprintf ( 1, '  SUM_SUB carries it out.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The integration interval is [%f, %f]\n', a, b );
  fprintf ( 1, '  Number of subintervals is %d\n', nsub );
  fprintf ( 1, '  Quadrature order will vary.\n' );
  fprintf ( 1, '  Integrand will vary.\n' );
  fprintf ( 1, '\n' );

  for ilo = 1 : 5 : nfunc

    ihi = min ( ilo + 4, nfunc );

    fprintf ( 1, '\n' );
    fprintf ( 1, '    ' );
    for i = ilo : ihi
      fprintf ( '%14s', fname(i) );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '\n' );

    for norder = 1 : order_max

      fprintf ( 1, '  %2d', norder );

      for i = ilo : ihi

        func_set ( 'SET', i );

        [ xtab, weight ] = legendre_set_x0_01 ( norder );

        result(i) = sum_sub ( @func, a, b, nsub, norder, xlo, xhi, ...
          xtab, weight );

        fprintf ( 1, '  %12f', result(i) );

      end

      fprintf ( 1, '\n' );

    end

  end

  return
end

author: johannesgerer


text:
function [y n] = var_nan(x)

for ii=1:size(x,2)
    tmp = x(:,ii);
    tmp(isnan(tmp))=[];
    y(ii)=var(tmp);
    n(ii)=size(tmp,1);
end


text:
function handles = togglePlotData(hObject, handles, state)

set(handles.btnPlotData, 'Enable', state);

if isequal(handles.foundExportFig, 1)
    set(handles.btnPlotDataSave, 'Enable', state);
end

guidata(hObject, handles);

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function varargout = myparse(params,varargin)

nin = length(varargin)/2;
nout = nargout;
if mod(2*nin,2)~=0 || mod(length(params),2)~=0,
  error('Odd number of parameters expected');
end;

% set defaults
varargout = varargin(2:2:2*nout);

% parameter names
names = varargin(1:2:end);

% go through parameters
for i = 1:2:length(params),

  % check if it matches names{j}
  if ~ischar(params{i}),
    error('Input %d is not a string.',i);
  end

  %BJA:  6x faster than that below
  idx=find(strcmpi(params{i},names));
  if(isempty(idx))
    warning('myparse:unknownParameterName','Unknown parameter name: %s, ignoring\n',params{i});
  else
    varargout{idx} = params{i+1};
  end

%  if ~any(strcmpi(params{i},names)),
%    warning('myparse:unknownParameterName','Unknown parameter name: %s, ignoring\n',params{i});
%  end
%
%  for j = 1:length(names),
%
%    % set value
%    if strcmpi(params{i},names{j}),
%      varargout{j} = params{i+1};
%    end;
%
%  end;

end;

author: kristinbranson

repo: JAABA


text:
function c = tapas_kf_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the Kalman filter
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The Kalman filter configuration consists of the priors of parameters and initial values. All
% priors are Gaussian in the space where the quantity they refer to is estimated. They are specified
% by their sufficient statistics: mean and variance (NOT standard deviation).
%
% Quantities are estimated in their native space if they are unbounded (e.g., mu_0). They are
% estimated in log-space if they have a natural lower bound at zero (e.g., the pi_u).
%
% Parameters can be fixed (i.e., set to a fixed value) by setting the variance of their prior to
% zero.
%
% Fitted trajectories can be plotted by using the command
%
% >> tapas_kf_plotTraj(est)
%
% where est is the stucture returned by tapas_fitModel. This structure contains the estimated
% perceptual parameters in est.p_prc and the estimated trajectories of the filter's
% representations. Their meanings are:
%
%         est.p_prc.g_0        initial value of gain
%         est.p_prc.mu_0       initial values of hidden state mean
%         est.p_prc.om         process variance
%         est.p_prc.pi_u       observation precision
%
%         est.traj.da          prediction error
%         est.traj.g           gain
%         est.traj.mu          hidden state mean
%
% Tips:
% - When analyzing a new dataset, take your inputs u and use
%
%   >> est = tapas_fitModel([], u, 'tapas_kf_config', 'tapas_bayes_optimal_config');
%
%   to determine the Bayes optimal perceptual parameters (given your current priors as defined in
%   this file here, so choose them wide and loose to let the inputs influence the result). You can
%   then use the optimal parameters as your new prior means for the perceptual parameters.
%
% - If the log-model evidence cannot be calculated because the Hessian poses problems, look at
%   est.optim.H and fix the parameters that lead to NaNs.
%
% - Your guide to all these adjustments is the log-model evidence (LME). Whenever the LME increases
%   by at least 3 across datasets, the adjustment was a good idea and can be justified by just this:
%   the LME increased, so you had a better model.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2016 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Config structure
c = struct;

% Model name
c.model = 'Kalman filter';

% Sufficient statistics of Gaussian parameter priors

% PLACEHOLDER VALUES
% It is often convenient to set some priors to values
% derived from the inputs. This can be achieved by
% using placeholder values. The available placeholders
% are:
%
% 99991   Value of the first input
%         Usually a good choice for mu_0mu(1)
% 99992   Variance of the first 20 inputs
%         Usually a good choice for mu_0sa(1)
% 99993   Log-variance of the first 20 inputs
%         Usually a good choice for logsa_0mu(1), and
%         its negative, ie the log-precision of the
%         first 20 inputs, for logpiumu
% 99994   Log-variance of the first 20 inputs minus two
%         Usually a good choice for ommu(1)

% Initial gain
c.logg_0mu = 0.1;
c.logg_0sa = 1;

% Initial hidden state mean
c.mu_0mu = 99991;
c.mu_0sa = 99992;

% Process variance
c.ommu = 99993;
c.omsa = 1;

% Pi_u
% Fix this to zero (no percpeptual uncertainty) by setting
% logpiumu = -Inf; logpiusa = 0;
c.logpiumu = -99993;
c.logpiusa = 1;

% Gather prior settings in vectors
c.priormus = [
    c.logg_0mu,...
    c.mu_0mu,...
    c.ommu,...
    c.logpiumu,...
         ];

c.priorsas = [
    c.logg_0sa,...
    c.mu_0sa,...
    c.omsa,...
    c.logpiusa,...
         ];

% Model function handle
c.prc_fun = @tapas_kf;

% Handle to function that transforms perceptual parameters to their native space
% from the space they are estimated in
c.transp_prc_fun = @tapas_kf_transp;

return;

author: translationalneuromodeling

repo: tapas


function [ber, rate] = OMP_method()

global Ns Nrf H Vn W_mopt Codebook_v Codebook_w;
n = 1;   %number of iterCodebook_vion
MSE =zeros(1,11);
W_D = W_mopt;  %initializCodebook_vion
W_RF = 1;
tw = trace(W_D'*(W_RF)'*W_RF*W_D);
[~,N_t] = size(H);
V_RF = exp( 1i*unifrnd(0,2*pi,N_t,Nrf) );
while(n<3 || (MSE(n-2)-MSE(n-1))>1e-4 &&n<=10)
    H_u = H'*W_RF*W_D;          %effective downlink channel
    Vn1 = tw*Vn;
    [V_RF,V_D] = MSEOMP (Nrf,H_u,Codebook_v,Vn1);
    tv = trace(V_RF*(V_D)*V_D'*V_RF');
    %%UEside
    H_d = H*V_RF*V_D;
    Vn2 = tv*Vn;
    [W_RF,W_D] = MSEOMP (Nrf,H_d,Codebook_w,Vn2); %the same formulCodebook_vion
    He = W_D'*W_RF'*H*V_RF*V_D;
    tw = trace(W_D'*(W_RF)'*W_RF*W_D);
    MSE(n) = trace(He * He'- He- He'+ eye(Ns)) + Vn2*tw;
    n = n + 1;
end
V_D = V_D/sqrt(tv);
%W_B = W_D*sqrt(tv);

V = V_RF * V_D;
W = W_RF * W_D;

ber = get_ber(V, W);
rate = get_rate(V, W);


text:
clear all; close all;
I=imread('peppers.png');
J=rgb2gray(I);
figure;
subplot(121);
imshow(J);
subplot(122);
imcontour(J,3);

author: HuangCongQing


text:
function r = exprand(mu,m,n)
% EXPRAND Random matrices from exponential distribution.
%
%    R = EXPRAND(MU) returns a matrix of random numbers chosen
%    from the exponential distribution with parameter MU.
%    The size of R is the size of MU.
%    Alternatively, R = EXPRAND(MU,M,N) returns an M by N matrix.

% Copyright (c) 1998-2004 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

error('No mex-file for this architecture. See Matlab help and convert.m in ./linuxCsource or ./winCsource for help.')

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [net_matching] = init_siamese_model()
% this file implements the matching network
% By Xin Li, 11-27-2017

net_matching = dagnn.DagNN();
    net_matching.addLayer('xcorr', XCorr(), ...
                 {'w_feat', 'x_feat'}, ...
                 {'m_score'}, ...
                 {});
net_matching.move('gpu');

end

author: XinLi-zn

repo: TADT

save_path: github-repos/MATLAB/XinLi-zn-TADT


text:
function yawm = mag_compass_update(m_n, dec, inc, DCMnb, roll, pitch)
% magh_update: calculates magnetic heading angle (yaw) from magnetometer data.
%
% INPUT
%   mag, 1x3 magnetic flux density (Tesla).
%   dec, magnetic declination angle (rad).
%   inc, magnetic inclination angle (rad).
%   DCMnb, 3x3 DCM nav-to-body.
%   roll, roll angle (rad).
%   pitch, pitch angle (rad).
%
% OUTPUT
%   yawm, yaw angle from magnetometer data.
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.
%
%   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 Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% Reference:
%
%   Paul D. Groves. Principles of GNSS, Inertial, and
% Multisensor Integrated Navigation Systems. Second Edition.
% Eq. 6.2 to 6.6, page 219.
%
%   NOAA, Magnetic Field Calculators.
% https://www.ngdc.noaa.gov/geomag/calculators/magcalc.shtml
%
% Version: 002
% Date:    2021/03/20
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

B = norm( m_n );                    % Magnitude of the flux density

D = [ cos(dec) * cos(inc); sin(dec) * cos(inc); sin(inc); ];

m_b = ( DCMnb * D * B ) ;

x =  -m_b(2) * cos(roll)  + m_b(3) * sin(roll) ;
y =   m_b(1) * cos(pitch) + m_b(2) * sin(roll) * sin(pitch) ...
    + m_b(3) * cos(roll) * sin(pitch) ;

yawm = correct_yaw (atan2(x , y) + dec);    % Four-quadrant arctangent function

end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
%% Low Level Parent Phase Reconstruction
%
%%
% In this section we discuss parent grain reconstruction at the example of
% a titanium alloy. Lets start by importing a sample data set

mtexdata alphaBetaTitanium

% the phase names for the alpha and beta phases
alphaName = 'Ti (alpha)';
betaName = 'Ti (Beta)';

% and plot the alpha phase as an inverse pole figure map
plot(ebsd(alphaName),ebsd(alphaName).orientations,'figSize','large')

%%
% The data set contains 99.8 percent alpha titanium and 0.2 percent beta
% titanium. Our goal is to reconstuct the original beta phase. The
% original grain structure appears almost visible for human eyes.
% Our computations will be based on the Burgers orientation relationship

beta2alpha = orientation.Burgers(ebsd(betaName).CS,ebsd(alphaName).CS)

%%
% that alligns (110) plane of the beta phase with the (0001) plane of the
% alpha phase and the [1-11] direction of the beta phase with the [2110]
% direction of the alpha phase.
%
% Note that all MTEX functions for parent grain reconstruction expect the
% orientation relationship as parent to child and not as child to parent.
%
%% Detecting triple points that belong to the same parent orientation
%
% In a first step we want to identify triple junctions that have
% misorientations that are compatible with a common parent orientations. To
% this end we first compute alpha grains using the option
% <QuadruplePoints.html |removeQuadruplePoints|> which turn all quadruple
% junctions into 2 triple junctions. Furthermore, we choose a very small
% threshold of 1.5 degree for the identification of grain boundaries to
% avoid alpha orientations that belong to different beta grains get merged
% into the same alpha grain.

% reconstruct grains
[grains,ebsd.grainId] = calcGrains(ebsd('indexed'),'threshold',1.5*degree,...
  'removeQuadruplePoints');
grains = smooth(grains,1,'moveTriplePoints');

% plot all alpha pixels
region = [299 401 -500 -440];
plot(ebsd(alphaName),ebsd(alphaName).orientations,...
  'region',region,'micronbar','off','figSize','large');

% and on top the grain boundaries
hold on
plot(grains.boundary,'linewidth',2 ,'region',region);
hold off

%%
% Above we have plotted only a very small subregion of the original data
% set to make the seperation of the qudruple junctions better visible.
%
% Next we extract all alpha - alpha - alpha triple junctions and use the
% command <calcParent.html |calcParent|> to find for each of these triple
% junctions the best fitting parent orientations.

% extract all alpha - alpha - alpha triple points
tP = grains.triplePoints(alphaName,alphaName,alphaName)

% compute for each triple point the best fitting parentId and how well the fit is
tPori = grains(tP.grainId).meanOrientation;
[parentId, fit] = calcParent(tPori,beta2alpha,'numFit',2,'id','threshold',5*degree);

%%
% The command |calcParent| returns for each child orientation a |parentId|
% which allows us later to compute the parent orientation from the child
% orientation. Furthermore, the command return for each triple junction the
% misfit between the adjecent parent orientations in radiant. Finally, the
% option |'numFit',2| causes |calcParent| to return not only the best fit
% but also the second best fit. This will be used later. First we simple
% colorize the triple junctions according to the best fit.

hold on
plot(tP,fit(:,1) ./ degree,'MarkerEdgecolor','k','MarkerSize',10,'region',region)
setColorRange([0,5])
mtexColorMap LaboTeX
mtexColorbar
hold off

%%
% Next we select those triple junctions as reliable that have a fit less
% than 2.5 degree and second best fit that is larger than  2.5 degree

consistenTP = fit(:,1) < 2.5*degree & fit(:,2) > 2.5*degree;

% marke these triple points by a red cicle
hold on
plot(tP(consistenTP),'MarkerEdgecolor','r','MarkerSize',10,...
  'MarkerFaceColor','none','linewidth',2,'region',region)
hold off

%% Recover beta grains from consistent triple junctions
%
% We observe that despite the quite sharp threshold we have many consistent
% triple points. In the next step we check wether all consistent triple
% junctions of a grain vote for the same parent orientation. Such a check
% for consistent votes can be computed by the command <majorityVote.html
% |majorityVote|> using the option |strict|.

% get a unique parentId vote for each grain
[parentId, numVotes] = majorityVote( tP(consistenTP).grainId, ...
  parentId(consistenTP,:,1), max(grains.id),'strict');

%%
% The command |majorityVote| returns for each grain with consistent
% parentId votes this unique parentId and for all other grains |NaN|.
% The second  output argument gives the number of these votes
%
% For all grains with at least 3 unique vote we now use the command
% <variants.html |variants|> to compute the parent orientation
% corresponding to the |parentId|. This parent orientations we assign as
% new |meanOrientation| to our grains.

% lets store the parent grains into a new variable
parentGrains = grains;

% change orientations of consistent grains from child to parent
parentGrains(numVotes>2).meanOrientation = ...
  variants(beta2alpha,grains(numVotes>2).meanOrientation,parentId(numVotes>2));

% update all grain properties that are related to the mean orientation
parentGrains = parentGrains.update;

%%
% Lets plot map of these reconstructed beta grains

% define a color key
ipfKey = ipfColorKey(ebsd(betaName));
ipfKey.inversePoleFigureDirection = vector3d.Y;

% and plot
plot(parentGrains(betaName), ...
  ipfKey.orientation2color(parentGrains(betaName).meanOrientation),'figSize','large')

%%
% We observe that this first step already results in very many Beta grains.
% However, the grain boundaries are still the boundaries of the original
% alpha grains. To overcome this, we merge all Beta grains that have a
% misorientation angle smaller then 2.5 degree.
%
% As an additional consistency check we verify that each parent
% grain has been reconstructed from at least 2 child grains. To this end we
% first make a testrun the merge operation and then revert all parent
% grains that that have less then two childs. This step may not nessesary
% in many case.

% test run of the merge operation
[~,parentId] = merge(parentGrains,'threshold',2.5*degree,'testRun');

% count the number of neighbouring child that would get merged with each child
counts = accumarray(parentId,1);

% revert all beta grains back to alpha grains if they would get merged with
% less then 1 other child grains
setBack = counts(parentId) < 2 & grains.phaseId == grains.name2id(alphaName);
parentGrains(setBack).meanOrientation = grains(setBack).meanOrientation;
parentGrains = parentGrains.update;

%%
% Now we perform the actual merge and the reconstruction of the parent
% grain boundaries.

% merge beta grains
[parentGrains,parentId] = merge(parentGrains,'threshold',2.5*degree);

% set up a EBSD map for the parent phase
parentEBSD = ebsd;

% and store there the grainIds of the parent grains
parentEBSD('indexed').grainId = parentId(ebsd('indexed').grainId);

plot(parentGrains(betaName), ...
  ipfKey.orientation2color(parentGrains(betaName).meanOrientation),'figSize','large')

%% Merge alpha grains to beta grains
%
% After the first two steps we have quite some alpha grains have not yet
% transformed into beta grains. In order to merge those left over alpha
% grains we check whether their misorientation with one of the neighbouring
% beta grains coincides with the parent to grain orientation relationship
% and if yes merge them evantually with the already reconstructed beta
% grains.
%
% First extract a list of all neighbouring alpha - beta grains

% all neighbouring alpha - beta grains
grainPairs = neighbors(parentGrains(alphaName), parentGrains(betaName));

%%
% and check how well they fit to a common parent orientation

% extract the corresponding meanorientations
oriAlpha = parentGrains( grainPairs(:,1) ).meanOrientation;
oriBeta = parentGrains( grainPairs(:,2) ).meanOrientation;

% compute for each alpha / beta pair of grains the best fitting parentId
[parentId, fit] = calcParent(oriAlpha,oriBeta,beta2alpha,'numFit',2,'id');

%%
% Similarly, as in the first step the command <calcParent.html
% |calcParent|> returns a list of |parentId| that allows the convert the
% child orientations into parent orientations using the command
% <variants.html |variants|> and the fitting to the given parent
% orientation. Similarly, as for the triple point we select only those
% alpha beta pairs such that the fit is below the threshold of 2.5 degree
% and at the same time the second best fit is above 2.5 degree.

% consistent pairs are those with a very small misfit
consistenPairs = fit(:,1) < 5*degree & fit(:,2) > 5*degree;

%%
% Next we compute for all alpha grains the majority vote of the surounding
% beta grains and change their orientation from alpha to beta

parentId = majorityVote( grainPairs(consistenPairs,1), ...
  parentId(consistenPairs,1), max(parentGrains.id));

% change grains from child to parent
hasVote = ~isnan(parentId);
parentGrains(hasVote).meanOrientation = ...
  variants(beta2alpha, parentGrains(hasVote).meanOrientation, parentId(hasVote));

% update grain boundaries
parentGrains = parentGrains.update;

% merge new beta grains into the old beta grains
[parentGrains,parentId] = merge(parentGrains,'threshold',5*degree);

% update grainId in the ebsd map
parentEBSD('indexed').grainId = parentId(parentEBSD('indexed').grainId);

% plot the result
color = ipfKey.orientation2color(parentGrains(betaName).meanOrientation);
plot(parentGrains(betaName),color,'linewidth',2)

%%
% The above step has merged

sum(hasVote)

%%
% alpha grains into the already reconstructed beta grain. This reduces the
% amount of grains not yet reconstructed to

sum(parentGrains('Ti (alpha').grainSize) ./ sum(parentGrains.grainSize)*100

%%
% percent. One way to proceed would be to repeat the steps of this section
% multiple time, maybe with increasing threshold, until the percentage of
% reconstructed beta grains is sufficiently high. Another approach in to
% consider the left over alpha grains as noise and use denoising techniques
% to replace them with beta orientations. This will be done in the last
% section.

%% Reconstruct beta orientations in EBSD map
%
% Until now we have only recovered the beta orientations as the mean
% orientations of the beta grains. In this section we want to set up the
% EBSD variable |parentEBSD| to contain for each pixel a reconstruction of
% the parent phase orientation.
%
% Therefore, we first identify all pixels that previously have been alpha
% titanium but now belong to a beta grain.

% consider only original alpha pixels that now belong to beta grains

isNowBeta = parentGrains.phaseId(max(1,parentEBSD.grainId)) == ebsd.name2id(betaName) &...
  parentEBSD.phaseId == ebsd.name2id(alphaName);

%%
% Next we can use once again the function <calcParent.html |calcParent|> to
% recover the original beta orientation from the measured alpha orientation
% giving the mean beta orientation of the grain.

% update beta orientation
[parentEBSD(isNowBeta).orientations, fit] = calcParent(parentEBSD(isNowBeta).orientations,...
  parentGrains(parentEBSD(isNowBeta).grainId).meanOrientation,beta2alpha);

%%
% We obtain even a measure |fit| for the corespondence between the beta
% orientation reconstructed for a single pixel and the beta orientation of
% the grain. Lets visualize this measure of fit

% the beta phase
plot(parentEBSD(isNowBeta),fit ./ degree,'figSize','large')
mtexColorbar
setColorRange([0,5])
mtexColorMap('LaboTeX')

hold on
plot(parentGrains.boundary,'lineWidth',2)
hold off

%%
% Lets finaly plot the reconstructed beta phase

plot(parentEBSD(betaName),ipfKey.orientation2color(parentEBSD(betaName).orientations),'figSize','large')

%% Denoising of the reconstructed beta phase
% As promised we end our discussion by applying denoising techniques to
% fill the remaining holes of alpha grains. To this end we first
% reconstruct grains from the parent orientations and throw away all small
% grains

[parentGrains,parentEBSD.grainId] = calcGrains(parentEBSD('indexed'),'angle',5*degree);

% remove all the small grains
parentEBSD = parentEBSD(parentGrains(parentGrains.grainSize > 15));

% redo grain reconstruction
[parentGrains,parentEBSD.grainId] = calcGrains(parentEBSD('indexed'),'angle',5*degree);

% smooth the grains a bit
parentGrains = smooth(parentGrains,5);

%%
% Finally, we denoise the remaining beta orientations and at the same time
% fill the empty holes. We choose a very small smoothing parameter |alpha|
% to keep as many details as possible.

F= halfQuadraticFilter;
F.alpha = 0.1;
parentEBSD = smooth(parentEBSD,F,'fill',parentGrains);

% plot the resulting beta phase
plot(parentEBSD(betaName),ipfKey.orientation2color(parentEBSD(betaName).orientations),'figSize','large')

hold on
plot(parentGrains.boundary,'lineWidth',3)
hold off

%%
% For comparison the map with original alpha phase and on top the recovered
% beta grain boundaries

plot(ebsd(alphaName),ebsd(alphaName).orientations,'figSize','large')

hold on
plot(parentGrains.boundary,'lineWidth',3)
hold off

%% Summary of relevant thresholds
%
% In parent grain reconstruction several parameters are involve are
% decicive for the success of the reconstruction
%
% * threshold for initial grain segmentation (1.5*degree)
% * maximum misfit at triple junctions (2.5 degree)
% * minimal misfit of the second best solution at triple junctions (2.5 degree)
% * minimum number of consistent votes (2)
% * threshold for merging beta grains (can be skipped)
% * threshold for merging alpha and beta grains (2.5 degree)
%%
%
% Visualize the misorientation to the mean reveals quite some fine
% structure in the reconstructed parent orientations.

cKey = axisAngleColorKey;
color = cKey.orientation2color(parentEBSD(betaName).orientations, parentGrains(parentEBSD(betaName).grainId).meanOrientation);
plot(parentEBSD(betaName),color)

hold on
plot(parentGrains.boundary,'lineWidth',3)
hold off

author: mtex-toolbox

repo: mtex

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


text:
function Z=setoperatorname(Z,name)

Z.extra.opname=name;


text:
function V = ReferenceVectorRegeneration(PopObj,V)
% Reference vector regeneration

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

    PopObj        = PopObj - repmat(min(PopObj,[],1),size(PopObj,1),1);
    [~,associate] = max(1-pdist2(PopObj,V,'cosine'),[],2);
    inValid       = setdiff(1:size(V,1),associate);
    V(inValid,:)  = rand(length(inValid),size(V,2)).*repmat(max(PopObj,[],1),length(inValid),1);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function y = eq(a, b, same)
%EQ Overloaded equality operator, or test for Layer instance equality
%   EQ(A, B), A == B returns a Layer that tests equality of the outputs of
%   two Layers (one of them may be constant).
%
%   EQ(A, B, 'sameInstance') checks if two variables refer to the same
%   Layer instance (i.e., calls the == operator for handle classes).

  if nargin <= 2
    y = Layer(@eq, a, b) ;
    y.numInputDer = 0 ;  % non-differentiable
  else
    assert(isequal(same, 'sameInstance'), ...
      'The only accepted extra flag for EQ is ''sameInstance''.') ;

    if ~isa(a, 'Layer') || ~isa(b, 'Layer')
      y = false ;
    else
      y = eq@handle(a, b) ;
    end
  end
end

repo: MFT


text:
%% gedge
% Below is a demonstration of the features of the |dihedralAngles| function

%%
clear; close all; clc;

%% Syntax
% |[h]=gedge(E,V,AE,5,AE);|

%% Description
% The |gedge| function aids in the visualization of edges e.g. of meshes.
% Edge coloring and transparency can be prescribed.
%
% SEe also |gpatch|.

%% Examples
%

%%
% Plot settings
fontSize=20;
edgeThickness=5;
cmap=gjet(250);

%%
% Create example data, mesh data for a hex mesh

%Create mesh for a cube
boxDim=2*ones(1,3);
boxEl=1*ones(1,3);
[meshStruct]=hexMeshBox(boxDim,boxEl);
Eh=meshStruct.E;
V=meshStruct.V;
F=meshStruct.F;

%Distorting the shape (shear)
d=eye(3,3); %Identity
d(1,2)=1; %Add shear
V=V*d; %Deform

%Get mesh edges
E=patchEdges(F,V);

%% Example 1: Visualizing edges in a single color

cFigure;
subplot(1,2,1); hold on;
title('gpatch visualization of mesh faces and edges');
gpatch(F,V,'rw','r',1,edgeThickness);
axisGeom; camlight headlight;

subplot(1,2,2); hold on;
title('gedge visualization of mesh edges only');
gedge(E,V,'r',edgeThickness);
axisGeom; camlight headlight;

drawnow;

%% Example 2: Using colormapping on the edges

%%
% Create example color data on edges, e.g. the mesh dihedral angle

[~,E,AE]=dihedralAngles(Eh,V,'hex8');
AE=180*(AE./pi);

%%

cFigure; hold on;
title('Color data on edges')
gedge(E,V,AE,edgeThickness);
axisGeom; camlight headlight;
colormap(cmap); colorbar;
gdrawnow;

%% Example 3: Specifying edge transparency as well

cFigure; hold on;
title('Alpha data and color data on edges')
gedge(E,V,AE,edgeThickness,AE);
axisGeom; camlight headlight;
colormap(cmap); colorbar;
alim([25 135]);
gdrawnow;

%% Example 4: Interpolated colors from vertex (node) data

[F,V]=geoSphere(1,1);
E=patchEdges(F);
VE=patchCentre(E,V);
CE=VE(:,1); %Example edge color data  = mid-edge x-coordinate
CV=V(:,1); %Example vertex color data = x-coordinate

%%

cFigure; hold on;
title('Color data on vertices')
h= gedge(E,V,CV,edgeThickness);
h.EdgeColor='interp'; %Turn on interpolated shading
axisGeom; camlight headlight;
colormap(cmap); colorbar;
gdrawnow;

%% Example 5: Specifying RGB color data on edges

C_rgb1=cmaperise(CE,gjet(250));
C_rgb2=cmaperise(CE,viridis(250));
C_rgb3=cmaperise(CE,spectral(250));

V2=V;
V2(:,1)=V2(:,1)+2.5;
V3=V2;
V3(:,1)=V3(:,1)+2.5;

%%

cFigure; hold on;
title('RGB data on edges')
h1= gedge(E,V ,C_rgb1,edgeThickness);
h2= gedge(E,V2,C_rgb2,edgeThickness);
h3= gedge(E,V3,C_rgb3,edgeThickness);
axisGeom; camlight headlight;
gdrawnow;

%% Example 6: Specifying RGB color data on vertices (nodes)

C_rgb1=cmaperise(CV,gjet(250));
C_rgb2=cmaperise(CV,viridis(250));
C_rgb3=cmaperise(CV,spectral(250));

%%

cFigure; hold on;
title('RGB data on vertices (interpolated)')
h1= gedge(E,V ,C_rgb1,edgeThickness); h1.EdgeColor='interp';
h2= gedge(E,V2,C_rgb2,edgeThickness); h2.EdgeColor='interp';
h3= gedge(E,V3,C_rgb3,edgeThickness); h3.EdgeColor='interp';
axisGeom; camlight headlight;
gdrawnow;

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

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_selectdata

% the issue is that ft_selectdata_new removes the dimord

data = [];
data.trial{1} = randn(2,100);
data.time{1}  = [0:99]./100;
data.label    = {'chan1';'chan2'};

data = ft_checkdata(data, 'datatype', 'timelock');

cfg = [];
cfg.channel = 'chan1';
data2 = ft_selectdata(cfg, data);

assert(isfield(data2,'dimord'));


text:
hydro = struct();

hydro = readWAMIT(hydro,'sphere.out',[]);
hydro = radiationIRF(hydro,15,[],[],[],[]);
hydro = radiationIRFSS(hydro,[],[]);
hydro = excitationIRF(hydro,15,[],[],[],[]);
writeBEMIOH5(hydro)
plotBEMIO(hydro)

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim