text:
function [] = anim8_DIC3DPP_faceMeasure_onImages_n_n(DIC3DPPresults,pairIndex,faceMeasureString,varargin)
%% function for plotting 3D-DIC post-processing results from step 4 projected on the 2D images
% plotting the images with the 3D face measure results plotted on top
% on the left side the images from the reference camera (reference image and current images), and on the right side the
% images from the deformed camera
%
% calling options:
% [] = anim8_DIC3D_faceMeasures_onImages_n_n(DIC3DPPresults,pairIndex,faceMeasureString);
% [] = anim8_DIC3D_faceMeasures_onImages_n_n(DIC3DPPresults,pairIndex,faceMeasureString,optStruct);
%
% INPUT:

%%
Points=DIC3DPPresults.DIC2Dinfo{pairIndex}.Points;
nCamRef=DIC3DPPresults.DIC2Dinfo{pairIndex}.nCamRef;
nCamDef=DIC3DPPresults.DIC2Dinfo{pairIndex}.nCamDef;
nImages=DIC3DPPresults.DIC2Dinfo{pairIndex}.nImages;
currentFacesLogic=DIC3DPPresults.FacePairInds==pairIndex;
currentPointIndex=find(DIC3DPPresults.PointPairInds==pairIndex);
firstCurrentPointIndex=currentPointIndex(1);
F=DIC3DPPresults.Faces(currentFacesLogic,:);
F=F-firstCurrentPointIndex+1;
FaceCorr=cell(nImages,1);
for ii=1:nImages
    FaceCorr{ii}=DIC3DPPresults.FaceCorrComb{ii}(currentFacesLogic,:);
end
ImSet=cell(2*nImages,1);
ImPaths=DIC3DPPresults.DIC2Dinfo{pairIndex}.ImPaths;

% if ImPaths are not valid (for example if using on another computer, ask
% user to provide a new folder where all the images are located.
try
    ImSet{1}=imread(ImPaths{1});
catch
    newPath=uigetdir([],'Path to images is invalid. Please provide the correct path to the images (the folder containing all camera folders with the processed gray images)');
    for ii=1:length(ImPaths)
        [a,b,c]=fileparts(ImPaths{ii});
        as = strsplit(a,{'\','/'});
        ImPaths{ii}=[newPath '\' as{end} '\' b c];
    end
end

for ii=1:2*nImages
    ImSet{ii}=imread(ImPaths{ii});
    if size(ImSet{ii},3)==3
        ImSet{ii}=rgb2gray(ImSet{ii});
    end
end

switch nargin
    case 3 % in case no results were entered
        optStruct=struct;
    case 4
        optStruct=varargin{1};
    otherwise
        error('wrong number of input arguments');
end

%% cut out point with large correlation coefficient
if ~isfield(optStruct,'CorCoeffCutOff')
    CorCoeffCutOff=max(max(cell2mat(FaceCorr)));
else
    CorCoeffCutOff=optStruct.CorCoeffCutOff;
end

for ii=1:nImages
    FaceCorr{ii}(FaceCorr{ii}>CorCoeffCutOff)=NaN;
end

%%
FC=cell(nImages,1);
switch faceMeasureString
    case {'J','Lamda1','Lamda2'}
        for ii=1:nImages
            FC{ii}=DIC3DPPresults.Deform.(faceMeasureString){ii}(currentFacesLogic,:);
        end
        if isfield(optStruct,'cMap')
            cMap=optStruct.cMap;
        else
            cMap=coldwarm;
        end
        if ~isfield(optStruct,'FClimits')
            FCmax=0;
            for ii=1:nImages
                FCmax=max([max(abs(FC{ii}(~isnan(FaceCorr{ii}))-1)) FCmax]);
            end
            FClimits=[1-FCmax 1+FCmax];
        else
            FClimits=optStruct.FClimits;
        end
    case {'Emgn','emgn','Eeq','eeq','EShearMax','eShearMax'}
        for ii=1:nImages
            FC{ii}=DIC3DPPresults.Deform.(faceMeasureString){ii}(currentFacesLogic,:);
        end
        if isfield(optStruct,'cMap')
            cMap=optStruct.cMap;
        else
            cMap='parula';
        end
        if ~isfield(optStruct,'FClimits')
            FCmax=0;
            for ii=1:nImages
                FCmax=max([max(FC{ii}(~isnan(FaceCorr{ii}))) FCmax]);
            end
            FClimits=[0 FCmax];
        else
            FClimits=optStruct.FClimits;
        end

    case {'Epc1','Epc2','epc1','epc2'}
        for ii=1:nImages
            FC{ii}=DIC3DPPresults.Deform.(faceMeasureString){ii}(currentFacesLogic,:);
        end
        if isfield(optStruct,'cMap')
            cMap=optStruct.cMap;
        else
            cMap=coldwarm;
        end
        if ~isfield(optStruct,'FClimits')
            FCmax=0;
            for ii=1:nImages
                FCmax=max([max(abs(FC{ii}(~isnan(FaceCorr{ii})))) FCmax]);
            end
            FClimits=[-FCmax FCmax];
        else
            FClimits=optStruct.FClimits;
        end
    otherwise
        error('unexpected face measure string. plots not created');
end

for ii=1:nImages
    FC{ii}(isnan(FaceCorr{ii}))=NaN;
end

%% cut out point with extreme face color values
if ~isfield(optStruct,'dataLimits')
    dataLimits=[min(min(cell2mat(FC))) max(max(cell2mat(FC)))];
else
    dataLimits=optStruct.dataLimits;
end

for ii=1:nImages
    FC{ii}(FC{ii}<dataLimits(1) | FC{ii}>dataLimits(2))=NaN;
end

%% plot
hf=cFigure;
hf.Units='normalized'; hf.OuterPosition=[.05 .05 .9 .9]; hf.Units='pixels';

% Ref
ii=1;
subplot(1,2,1)
hp1=imagesc(repmat(ImSet{ii},1,1,3)); hold on;
hp2=gpatch(F,Points{ii},FC{ii},'none',0.5);
pbaspect([size(ImSet{ii},2) size(ImSet{ii},1) 1])
hs1=title(['Ref (Cam ' num2str(nCamRef) ' frame ' num2str(1) ')']);
colormap(cMap);
hc1=colorbar;
caxis(FClimits)
title(hc1, faceMeasureString);
hc1.FontSize=16;
axis off

% Cur
ii=nImages+1;
subplot(1,2,2)
hp3=imagesc(repmat(ImSet{ii},1,1,3)); hold on
hp4=gpatch(F,Points{ii},FC{1},'none',0.5);
pbaspect([size(ImSet{ii},2) size(ImSet{ii},1) 1])
hs2=title(['Cur ' num2str(ii) ' (Cam ' num2str(nCamDef) ' frame ' num2str(1) ')']);
colormap(cMap);
hc2=colorbar;
caxis(FClimits);
title(hc2, faceMeasureString);
hc2.FontSize=16;
axis off

drawnow

%Create the time vector
animStruct.Time=linspace(0,1,nImages);

for ii=1:nImages
    Pnow1=Points{ii};
    Pnow2=Points{ii+nImages};

    cNow1=FC{ii};
    cNow2=FC{ii};

    TitleNow1=['Cur ' num2str(ii) ' (Cam ' num2str(nCamRef) ' frame ' num2str(ii) ')'];
    TitleNow2=['Cur ' num2str(ii) ' (Cam ' num2str(nCamDef) ' frame ' num2str(ii) ')'];

    %Set entries in animation structure
    animStruct.Handles{ii}=[hp1,hp3,hp2,hp2,hp4,hp4,hs1,hs2]; %Handles of objects to animate
    animStruct.Props{ii}={'CData','CData','Vertices','CData','Vertices','CData','String','String'}; %Properties of objects to animate
    animStruct.Set{ii}={repmat(ImSet{ii},1,1,3),repmat(ImSet{ii+nImages},1,1,3),Pnow1,cNow1,Pnow2,cNow2,TitleNow1,TitleNow2}; %Property values for to set in order to animate

end

anim8(hf,animStruct);
addFigureButtons;
end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function [regions] = selectRegions(image_data)

BND_THRESH = 0.010;
AREA_THRESH = 20*20;

pb1 = image_data.occ.pb1;
pb2 = image_data.occ.pb2;
bndinfo = image_data.occ.bndinfo_all{1};

hier = boundaries2hierarchy(pb1+pb2, bndinfo.edges.spLR, 'mean');
cost = [hier.init_cost ; hier.cost];
stats = regionprops(bndinfo.wseg, 'Area');
area = cat(1, stats.Area);
keep = false(size(hier.regions));

for r = 1:numel(keep)
   keep(r) = ((sum(area(hier.regions{r}))>AREA_THRESH) && (cost(r)>BND_THRESH));
end
regions = hier.regions(keep);

author: Cloud-CV

repo: object-proposals


text:
function [Bdraw,statOK]=sampleB(yData,Psi,SV_H,A,p,priorValues)
% This function uses the method proposed by Carriero Clark and Marcellino
% to draw from the conditional posterior distribution of B.
%sampleB(yData,PsiDraw_prop,HvarsDraw_prop,Adraw_prop,p,priorValues)
%Psi = PsiDraw_prop; %local mean
%SV_H = HvarsDraw_prop; %draw for the time varrying diagonal elements of the VCV
%A = Adraw_prop;     %the time invariant component of the VCV

%% Initialize
[T,M] = size(yData);
Tp=T-p;

% expand prior data
varsH=priorValues.vars;
lambda=priorValues.lambda; %overall tightness
theta=priorValues.theta;   %extra shrinkage for off diagonal elements

%% set up prior variance and mean of B
% Prior mean of B
priorMeanBeta=zeros(M*p,1); %means are equal to zero

% Prior variance of B
priorVarVectors=zeros(M*p,p);

% lag shrinkage
pVector=zeros(p,1);
for j=1:p
    pVector(j)=1/(j^2);
end

varsH=reshape(varsH,M,1);

for i=1:M
    % scale by variance estimates
    vector_iTemp=(varsH(i)*ones(M,1))./varsH;
    vector_iTemp=vector_iTemp*lambda*theta;
    vector_iTemp(i)=vector_iTemp(i)/theta;

    % calculate the prior variances
    vector_i=kron(pVector,vector_iTemp);
    priorVarVectors(:,i)=vector_i;
end

%% Now sample B row-by-row
Ainv=A\eye(M);

Y_Psi=yData-Psi;
Y_Psi(1:p,:)=yData(1:p,:)-ones(p,1)*mean(yData(1:p,:));
X_Psi = lagmatrix(Y_Psi,1:p);
X_Psi = X_Psi(p+1:end,:);
Y_Psi=Y_Psi(p+1:end,:);

OK=0; % to control for unstable draws

rejections=0;
statOK=1;
% start drawing B
while OK==0

    % initialize
    scaledEpsilon=zeros(Tp,M);
    Bdraw=zeros(p*M,M);

    % prepare data for column 1
    y1=Y_Psi(:,1);
    y1Scaled=y1./(SV_H(p+1:T,1).^0.5);
    hScaling1=(SV_H(p+1:T,1).^0.5)*ones(1,p*M);
    X1Scaled=X_Psi./hScaling1;

    % calculate posterior distribution
    postVarBinv=diag(priorVarVectors(:,1).^(-1))+X1Scaled'*X1Scaled;
    postVarB=postVarBinv\eye(M*p);
    postMeanB=postVarB*(diag(priorVarVectors(:,1).^(-1))*priorMeanBeta+X1Scaled'*y1Scaled);

    % obtain cholvar
    [cholPostVarBi,testPD]=chol(postVarB);
    if testPD>0
        cholPostVarBi= cholred(postVarB);
        disp('NPD!')
    end
    cholPostVarBi=cholPostVarBi'; % transpose to obtain lower triangular matrix

    % sample the column of Bdraw
    Bi_draw=postMeanB+cholPostVarBi*randn(p*M,1);
    Bdraw(:,1)=Bi_draw;

    % prepare residuals
    resids_1=y1Scaled-X1Scaled*Bi_draw;
    scaledEpsilon(:,1)=resids_1;

    for i=2:M

        % prepare data for column i
        a_vector=Ainv(i,:)';
        yi=Y_Psi(:,i)-scaledEpsilon*a_vector;
        yiScaled=yi./(SV_H(p+1:T,i).^0.5);

        hScaling_i=(SV_H(p+1:T,i).^0.5)*ones(1,p*M);
        XiScaled=X_Psi./hScaling_i;

        % calculate posterior distribution
        postVarBinv=diag(priorVarVectors(:,i).^(-1))+XiScaled'*XiScaled;
        postVarB=postVarBinv\eye(M*p);
        postMeanB=postVarB*(diag(priorVarVectors(:,i).^(-1))*priorMeanBeta+XiScaled'*yiScaled);

        % obtain cholvar
        [cholPostVarBi,testPD]=chol(postVarB);
        if testPD>0
            cholPostVarBi= cholred(postVarB);
            disp('NPD!')
        end
        cholPostVarBi=cholPostVarBi'; % transpose to obtain lower triangular matrix

        % sample the column of Bdraw
        Bi_draw=postMeanB+cholPostVarBi*randn(p*M,1);
        Bdraw(:,i)=Bi_draw;

        % prepare residuals
        resids_i=yiScaled-XiScaled*Bi_draw;
        scaledEpsilon(:,i)=resids_i;
    end

    % check stability
    [ max_vTemp ] = determineEV(Bdraw,M,p);
    max_v=abs(max_vTemp);

    if max_v <0.999
        OK=1; % accept the draw
    else
        rejections=rejections+1;
        %disp([max_v rejections])
        if rejections>1000
%             disp('Redo This round')
            statOK=0;
            OK=1;
        end

    end

end

end

author: european-central-bank

repo: BEAR-toolbox


text:
function varargout = pr(varargin)
%VL_PR   Precision-recall curve.
%   [RECALL, PRECISION] = VL_PR(LABELS, SCORES) computes the
%   precision-recall (PR) curve. LABELS are the ground truth labels,
%   greather than zero for a positive sample and smaller than zero for
%   a negative one. SCORES are the scores of the samples obtained from
%   a classifier, where lager scores should correspond to positive
%   samples.
%
%   Samples are ranked by decreasing scores, starting from rank 1.
%   PRECISION(K) and RECALL(K) are the precison and recall when
%   samples of rank smaller or equal to K-1 are predicted to be
%   positive and the remaining to be negative. So for example
%   PRECISION(3) is the percentage of positive samples among the two
%   samples with largest score. PRECISION(1) is the precision when no
%   samples are predicted to be positive and is conventionally set to
%   the value 1.
%
%   Set to zero the lables of samples that should be ignored in the
%   evaluation. Set to -INF the scores of samples which are not
%   retrieved. If there are samples with -INF score, then the PR curve
%   may have maximum recall smaller than 1, unless the INCLUDEINF
%   option is used (see below). The options NUMNEGATIVES and
%   NUMPOSITIVES can be used to add additional surrogate samples with
%   -INF score (see below).
%
%   [RECALL, PRECISION, INFO] = VL_PR(...) returns an additional
%   structure INFO with the following fields:
%
%   info.auc::
%     The area under the precision-recall curve. If the INTERPOLATE
%     option is set to FALSE, then trapezoidal interpolation is used
%     to integrate the PR curve. If the INTERPOLATE option is set to
%     TRUE, then the curve is piecewise constant and no other
%     approximation is introduced in the calculation of the area. In
%     the latter case, INFO.AUC is the same as INFO.AP.
%
%   info.ap::
%     Average precision as defined by TREC. This is the average of the
%     precision observed each time a new positive sample is
%     recalled. In this calculation, any sample with -INF score
%     (unless INCLUDEINF is used) and any additional positive induced
%     by NUMPOSITIVES has precision equal to zero. If the INTERPOLATE
%     option is set to true, the AP is computed from the interpolated
%     precision and the result is the same as INFO.AUC. Note that AP
%     as defined by TREC normally does not use interpolation [1].
%
%   info.ap_interp_11::
%     11-points interpolated average precision as defined by TREC.
%     This is the average of the maximum precision for recall levels
%     greather than 0.0, 0.1, 0.2, ..., 1.0. This measure was used in
%     the PASCAL VOC challenge up to the 2008 edition.
%
%   info.auc_pa08::
%     Deprecated. It is the same of INFO.AP_INTERP_11.
%
%   VL_PR(...) with no output arguments plots the PR curve in the
%   current axis.
%
%   VL_PR() accepts the following options:
%
%   Interpolate:: false
%     If set to true, use interpolated precision. The interpolated
%     precision is defined as the maximum precision for a given recall
%     level and onwards. Here it is implemented as the culumative
%     maximum from low to high scores of the precision.
%
%   NumPositives:: []
%   NumNegatives:: []
%     If set to a number, pretend that LABELS contains this may
%     positive/negative labels. NUMPOSITIVES/NUMNEGATIVES cannot be
%     smaller than the actual number of positive/negative entrires in
%     LABELS. The additional positive/negative labels are appended to
%     the end of the sequence, as if they had -INF scores (not
%     retrieved). This is useful to evaluate large retrieval systems
%     for which one stores ony a handful of top results for efficiency
%     reasons.
%
%   IncludeInf:: false
%     If set to true, data with -INF score SCORES is included in the
%     evaluation and the maximum recall is 1 even if -INF scores are
%     present. This option does not include any additional positive or
%     negative data introduced by specifying NUMPOSITIVES and
%     NUMNEGATIVES.
%
%   Stable:: false
%     If set to true, RECALL and PRECISION are returned in the same order
%     of LABELS and SCORES rather than being sorted by decreasing
%     score (increasing recall). Samples with -INF scores are assigned
%     RECALL and PRECISION equal to NaN.
%
%   NormalizePrior:: []
%     If set to a scalar, reweights positive and negative labels so
%     that the fraction of positive ones is equal to the specified
%     value. This computes the normalised PR curves of [2]
%
%   About the PR curve::
%     This section uses the same symbols used in the documentation of
%     the VL_ROC() function. In addition to those quantities, define:
%
%       PRECISION(S) = TP(S) / (TP(S) + FP(S))
%       RECALL(S) = TPR(S) = TP(S) / P
%
%     The precision is the fraction of positivie predictions which are
%     correct, and the recall is the fraction of positive labels that
%     have been correctly classified (recalled). Notice that the recall
%     is also equal to the true positive rate for the ROC curve (see
%     VL_ROC()).
%
%   REFERENCES:
%   [1] C. D. Manning, P. Raghavan, and H. Schutze. An Introduction to
%   Information Retrieval. Cambridge University Press, 2008.
%   [2] D. Hoiem, Y. Chodpathumwan, and Q. Dai. Diagnosing error in
%   object detectors. In Proc. ECCV, 2012.
%
%   See also VL_ROC(), VL_HELP().
[varargout{1:nargout}] = vl_pr(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
%  Figure 10.42      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
% fig10_42.m is a script to generate Figure 10.42, the root locus for
% altitude control with inner-loop stabilization and altitude feedback
% alone

clf;
ftq = [-0.0064    0.0263         0  -32.2000         0;
   -0.0941   -0.6240  761.1400 -196.2000         0;
   -0.0002   -0.0015   -4.4120  -12.4800         0;
         0         0    1.0000         0         0;
         0   -1.0000         0  830.0000         0];

g =[0;
  -32.7000;
   -2.0800;
         0;
         0];

h=[0 0 0 0 1];

j=0;

rlocus(ftq,g,-h,j);
v=[-12 12 -4 4];
axis(v);
grid;
title('Fig. 10.42  Rootlocus for altitude feedback with inner-loop stab.')

author: Sable

repo: mcbench-benchmarks


text:
function [data,units] = compute_max_wing_angle(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  % signed maximum: left wing will be negative, right wing positive
  data{i} = max(-trx(fly).wing_anglel,trx(fly).wing_angler);

end
units = parseunits('rad');

author: kristinbranson

repo: JAABA


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: computes the angular displacements of a damped pendulum under
%           gravity
%
% PURPOSE: to illustrate that changing the damping parameter, b,
%          significantly affects the resulting dynamics of the pendulum system
%
% Author: N.A. Battista
% Date: 12/13/2019
% Institution: The College of New Jersey
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function damped_pendulum_diff_b()

%
% Physical Parameters
%
m = 100;    % Mass of Pendulum bob
g=9.80665;  % Gravitational Acceleration;
L=0.2;      % Length of Pendulum

%
% Radii of Bobs (here - all the same)
%
rVec = 100*0.005*ones(1,6);

%
% Damping Parameter (here - each case is different)
%
bVec = 100*[0 0.1 0.2 0.4 1.0 2.0];

%
% Temporal Parameters
%
t=0;        % initial time
tFinal = 8; % final-time
dt = 1e-3;  % time-step
ct=1;       % counter
tVec(ct)=t; % time vector storage

%
% Moment of Inertia
%
IVec = m*(0.5*rVec.^2+L^2);

%
% Initial Angular Displacement
%
ang = -(pi/2-pi/5)*ones(1,6);
phi = zeros( size( ang ) );

%
% Perform the time-stepping, e.g., solve the ODEs
%
while t<=tFinal

    % increment the storage counter
    ct = ct+1;

    % increment time
    t = t+dt;

    % save current time
    tVec(ct) = t;

    % solve 1st order ODE for angular displacement
    ang(ct,:) = ang(ct-1,:) + dt * ( phi(ct-1,:) );

    % solve 1st order ODE for angular velocity
    phi(ct,:) = phi(ct-1,:) + dt * ( -m*g*L ./IVec .* sin( ang(ct-1,:) ) - bVec.*phi(ct-1,:)./IVec );

end

%
% Defining Colors for Plotting Data
%
color2 = [0.635 0.078 0.184];  % dark red
color5 = [0.9 0.525 0.098];    % orange
color6 = [0.929 0.840 0.1250]; % dark yellow
color7 = [0.7 0.7 0];          % 'ugly' green
color13 = [0 0.5 0.5];             % green

%
% Note discrepancy is a bit off due to small angle approximation:
% sin(theta) ~ theta gives below period of oscillation
%
T1 = 2*pi*sqrt( IVec(1) / (m*g*L) );
%
% Actually computed value of period is:
T1 = 1.93;

%
% Make Figure of Angular Displacement vs. Non-dimensional Time (# of oscillations of 'r' case)
%
lw=6;
fs=18;
plot(tVec/T1,ang(:,1),'k-','LineWidth',lw); hold on;
plot(tVec/T1,ang(:,2),'b-','LineWidth',lw,'Color',color2); hold on;
plot(tVec/T1,ang(:,3),'g-','LineWidth',lw,'Color',color5); hold on;
plot(tVec/T1,ang(:,4),'r-','LineWidth',lw,'Color',color6); hold on;
plot(tVec/T1,ang(:,5),'c-','LineWidth',lw,'Color',color7); hold on;
plot(tVec/T1,ang(:,6),'c-','LineWidth',lw,'Color',color13); hold on;
xlabel('Non-Dimensional Time (# periods of r-case)');
ylabel('Angular Displacement (Radians)');
axis([0 4.1 -1 1.5]);
leg=legend('b=0','b=10','b=20','b=40','b=100','b=200','Orientation','horizontal');
set(gca,'FontSize',fs);
set(leg,'FontSize',fs);
leg.NumColumns=3;

author: nickabattista


text:
%[ P, inls ] = ht_simple_ransac_p3p( u, X, rthr, maxiter )
%u: 3 x n image points
%X: 3 x n 3D points
%rthr: inlier threshold
%maxiter: default=1000

function [ P, inls ] = ht_simple_ransac_p3p( u, X, rthr, max_iter )
if nargin < 4
    max_iter = 1000;
end

%initialization
u = bsxfun(@rdivide, u, sqrt(sum(u.^2, 1)));
Npts = size(u, 2);
rthr = cos(rthr);
max_inlsnum = 3;
no_iter = 0;
P = [];
inls = false(1, Npts);
%ransac
while no_iter < max_iter
    no_iter = no_iter + 1;

    idx = randperm(Npts, 3);
    P_cand = P3PSolver([u(:, idx); X(:, idx)]);
    [inls_cand, inls_cand_num] = calculate_inls_angular(P_cand, u, X, rthr);
    if length(P_cand) > 1
        [inls_cand_num, inls_cand_idx] = max(inls_cand_num);
        inls_cand = inls_cand{inls_cand_idx};
        P_cand = P_cand{inls_cand_idx};
    else
        inls_cand_num = inls_cand_num(1);
        inls_cand = inls_cand{1};
        P_cand = P_cand{1};
    end

    if inls_cand_num >= max_inlsnum
        max_inlsnum = inls_cand_num;
        P = P_cand;
        inls = inls_cand;
        max_iter = min([max_iter, nsamples(max_inlsnum, Npts, 3, 0.95)]);
    end

end

end

function [SampleCnt, q] = nsamples(ni, ptNum, pf, conf)
    q  = prod (((ni-pf+1) : ni) ./ ((ptNum-pf+1) : ptNum));
    if q < eps
      SampleCnt = Inf;
    else
      %   SampleCnt  = log(1 - conf) / log(1 - q);
      if q > conf
        SampleCnt = 1;
      else
        SampleCnt  = log(1 - conf) / log(1 - q);
      end
    end
end

function [inls, inls_num] = calculate_inls_angular(Pcand, u, X, rthr)
    inls = cell(1, length(Pcand));
    inls_num = zeros(1, length(Pcand));
    for ii = 1:1:length(Pcand)
        X_reproj = Pcand{ii} * [X; ones(1, size(X, 2))];
        X_reproj = bsxfun(@rdivide, X_reproj, sqrt(sum(X_reproj.^2, 1)));

        res = sum(u .* X_reproj, 1);
        inls{ii} = res > rthr;
        inls_num(ii) = sum(inls{ii});
    end

end

author: HajimeTaira


text:
% Takes a hsvargplvm model and adds a non-(standard normal) prior on the
% parent (see: addDynamics functions in the other vargplvm-related
% packages).

% See also: svargplvmAddDynamics.m

function model = hsvargplvmAddParentPrior(model, globalOpt, optionsDyn)

modelParent = model.layer{model.H};
modelParent.comp{1}.vardist = modelParent.vardist;
modelParent.comp{1}.X = modelParent.vardist.means;
% TODO: if center means....
if model.H > 1
    modelParent.comp{1}.y = model.layer{end-1}.vardist.means;
elseif length(model.layer{1}.comp) == 1
    modelParent.comp{1}.y = model.layer{1}.comp{1}.y;
else
    error('Not implemented!')
end

modelParent.comp{1} = svargplvmAddDynamics(modelParent.comp{1}, globalOpt, optionsDyn);

modelParent.vardist = modelParent.comp{1}.vardist;
modelParent.dynamics = modelParent.comp{1}.dynamics;
modelParent.comp{1} = rmfield(modelParent.comp{1}, 'vardist');
modelParent.comp{1} = rmfield(modelParent.comp{1}, 'dynamics');
modelParent.comp{1} = rmfield(modelParent.comp{1}, 'X');
modelParent.comp{1} = rmfield(modelParent.comp{1}, 'y');

model.layer{model.H} = modelParent;

%params = hsvargplvmExtractParam(model);
%model = hsvargplvmExpandParam(model, params);

author: SheffieldML

repo: deepGP

save_path: github-repos/MATLAB/SheffieldML-deepGP


text:
%IPAD Pad an image with constants
%
% OUT = IPAD(IM, SIDES, N) is a padded version of the image IM with a block
% of NaN values N pixels wide on the sides of IM as specified by SIDES.
%
% OUT = IPAD(IM, SIDES, N, V) as above but pads with pixels of value V.
%
% SIDES is a string containing one or more of the characters:
% 't'   top
% 'b'   bottom
% 'l'   left
% 'r'   right
%
% Examples::
%
% Add a band of zero pixels 20 pixels high across the top of the image:
%     ipad(im, 't', 20, 0)
%
% Add a band of white pixels 10 pixels wide on all sides of the image:
%     ipad(im, 'tblr', 10, 255)
%
% Notes::
% - Not a tablet computer.

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

function out = ipad(in, sides, n, val)

    if nargin < 4
        val = NaN;
    end

    if ndims(in) > 2
        % multiplane case
        d = size(in);
        nimages = prod(d(3:end));
        out = [];
        for i=1:nimages
            out(:,:,i) = ipad(in(:,:,i), sides, n, val);
        end

        d2 = size(out);
        out = reshape(out, [d2(1) d2(2) d(3:end)]);
        return
    end

    out = in;
    for side=sides
        [w,h] = isize(out);

        switch side
            case 't'
                out = [val*ones(n,w); out];
            case 'b'
                out = [out; val*ones(n,w)];
            case 'l'
                out = [val*ones(h,n) out];
            case 'r'
                out = [out val*ones(h,n)];
        end
    end

author: petercorke

repo: machinevision-toolbox-matlab


text:
% -------------------------------------------------------------------------- %
% function hfssPolyline(fid, Name, Points, Units, [Closed], [segmentType], ...
%                       [Color], [Transparency])
% Description:
% ------------
%
% Parameters:
% -----------
% Name   - Name Attribute for the PolyLine.
% Points - Points as 3-Tuples, ex: Points = [0, 0, 1; 0, 1, 0; 1, 0 1];
%          Note: size(Points) must give [nPoints, 3].
% Units  - units of the segments (specify using either 'in', 'mm', 'meter'
%          or anything else defined in HFSS).
% Closed  - (Optional; bool) Closes the polyline to form a sheet. By
%           default is false.
% segmentType - (Optional) Defines the cross-section of the segment. Only
%               'Line' is supported. By default is 'Line'.
% Color       - (Optional) Defines the RGB color with a 3-tuple [R, G, B].
%               By default is set to [0, 0, 0] (black).
% Transparency - (Optional) Defines the transparency of the line. By
%                default is set to 0.8.
%
% Example:
% --------
% -------------------------------------------------------------------------- %

% ----------------------------------------------------------------------------
% This file is part of HFSS-MATLAB-API.
%
% HFSS-MATLAB-API 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.
%
% HFSS-MATLAB-API 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
% Foobar; if not, write to the Free Software Foundation, Inc., 59 Temple
% Place, Suite 330, Boston, MA  02111-1307  USA
%
% Copyright 2004, Vijay Ramasami (rvc@ku.edu)
% ----------------------------------------------------------------------------

% ----------------------------------------------------------------------------
% CHANGELOG
%
% ??-???-????: *Initial release (VR).
% 07-Aug-2014: *Fix representation of polylines (DRP).
% 07-Aug-2014: *Added option for closed polyline (DRP).
% 01-Sep-2020: *Fixed closed polyline assignment (DRP).
% 20-Apr-2021: *Added missing documentation (DRP).
% 21-Apr-2021: *Changed default color and transparency to that of the
%               hfssRectangle function (DRP).
% ----------------------------------------------------------------------------

% ----------------------------------------------------------------------------
% Modified by Daniel Rodriguez Prado
% danysan@gmail.com
% 20 Apr 2021
% ----------------------------------------------------------------------------

function hfssPolyline(fid, Name, Points, Units, Closed, segmentType, ...
                      Color, Transparency)
if (nargin < 5)
    Closed = [];
	segmentType = [];
	Color = [];
	Transparency = [];
elseif (nargin < 6)
    segmentType = [];
	Color = [];
	Transparency = [];
elseif (nargin < 7)
    Color = [];
	Transparency = [];
elseif (nargin < 8)
    Transparency = [];
end

if isempty(Closed)
    Closed = 'false';
end
if isempty(segmentType)
	segmentType = 'Line';
end
if isempty(Color)
	Color = [132, 132, 193];
end
if isempty(Transparency)
	Transparency = 0.75;
end

nPoints = length(Points(:, 1));

% Preamble.
fprintf(fid, '\n');
fprintf(fid, 'oEditor.CreatePolyline _\n');
fprintf(fid, '\tArray("NAME:PolylineParameters", ');
fprintf(fid, '"IsPolylineCovered:=", true, ');
if (Closed == true)
    fprintf(fid, '"IsPolylineClosed:=", true, _\n');
else
    fprintf(fid, '"IsPolylineClosed:=", false, _\n');
end

% Enter the Points.
fprintf(fid, '\t\tArray("NAME:PolylinePoints", _\n');
for iP = 1:nPoints-1
	fprintf(fid, '\t\t\tArray("NAME:PLPoint", ');
	fprintf(fid, '"X:=", "%.4f%s", ', Points(iP, 1), Units);
	fprintf(fid, '"Y:=", "%.4f%s", ', Points(iP, 2), Units);
	fprintf(fid, '"Z:=", "%.4f%s"), _\n', Points(iP, 3), Units);
end
fprintf(fid, '\t\t\tArray("NAME:PLPoint", ');
fprintf(fid, '"X:=", "%.4f%s",   ', Points(nPoints, 1), Units);
fprintf(fid, '"Y:=", "%.4f%s",   ', Points(nPoints, 2), Units);
fprintf(fid, '"Z:=", "%.4f%s")_\n', Points(nPoints, 3), Units);
fprintf(fid, '\t\t\t), _ \n');

% Create Segments.
fprintf(fid, '\t\tArray("NAME:PolylineSegments", _\n');
for iP = 1:nPoints-2
    fprintf(fid, '\t\t\tArray("NAME:PLSegment", ');
    fprintf(fid, '"SegmentType:=",  "%s", ', segmentType);
    fprintf(fid, '"StartIndex:=", %d, ', iP-1);
    fprintf(fid, '"NoOfPoints:=", 2), _\n');
end
fprintf(fid, '\t\t\tArray("NAME:PLSegment", ');
fprintf(fid, '"SegmentType:=",  "%s", ', segmentType);
fprintf(fid, '"StartIndex:=", %d, ', iP);
fprintf(fid, '"NoOfPoints:=", 2) _\n');
fprintf(fid, '\t\t\t) _\n');
fprintf(fid, '\t\t), _\n');

% Polyline Attributes.
fprintf(fid, 'Array("NAME:Attributes", _\n');
fprintf(fid, '"Name:=", "%s", _\n', Name);
fprintf(fid, '"Flags:=", "", _\n');
fprintf(fid, '"Color:=", "(%d %d %d)", _\n', Color(1), Color(2), Color(3));
fprintf(fid, '"Transparency:=", %f, _\n', Transparency);
fprintf(fid, '"PartCoordinateSystem:=", "Global", _\n');
fprintf(fid, '"MaterialName:=", "vacuum", _\n');
fprintf(fid, '"SolveInside:=", true)\n');

author: yuip

repo: hfss-api

save_path: github-repos/MATLAB/yuip-hfss-api


text:
function [uOutput] = cnssimp(nFunction, sFilename)
    % CNSSIMP import data of the ANSS/CNSS raw format catalog
    %
% [uOutput] = cnssimp(nFunction, sFilename);
%
% ANSS/CNSS raw format: (http://quake.geo.berkeley.edu/ncedc/documents.html#catalog_formats)
% Description of the parameters can be found at http://quake.geo.berkeley.edu/ftp/pub/doc/cat5/cnss.catalog.5
%
% updated: 15.09.03, J. Woessner

if nFunction == FilterOp.getDescription
    uOutput = 'ANSS/CNSS format (string conversion)';
elseif nFunction == FilterOp.importCatalog

    mData = textread(sFilename, '%s', 'delimiter', '\n', 'whitespace', '');

    %transform data to ZMAP format
    uOutput = zeros(size(mData, 1), 9);

    for i = 1:length(mData)
        if rem(i,100) == 0 ;
            disp([ num2str(i) ' of ' num2str(length(mData)) ' events processed ']);
        end
        try
            %             uOutput(i,1) = str2num(mData{i}(33:41));    %lon
            %             uOutput(i,2) = str2num(mData{i}(24:31));    %lat
            %             uOutput(i,3) = str2num(mData{i}(1:4));      %yr
            %             uOutput(i,4) = str2num(mData{i}(6:7));      %mo
            %             uOutput(i,5) = str2num(mData{i}(9:10));     %da
            %             uOutput(i,6) = str2num(mData{i}(51:54));    %mag

            uOutput(i,1) = str2num(mData{i}(34:43));    %lon
            uOutput(i,2) = str2num(mData{i}(25:33));    %lat
            uOutput(i,3) = str2num(mData{i}(6:9));      %yr
            uOutput(i,4) = str2num(mData{i}(10:11));      %mo
            uOutput(i,5) = str2num(mData{i}(12:13));     %da
            uOutput(i,6) = str2num(mData{i}(130:134));    %mag
            str = '      ';
            if  strcmp(mData{i}(44:51),str)== 1%strcmp(mData{i}(43:48),str)== 1
                uOutput(i,7) = 0;
            else
                %uOutput(i,7) = str2num(mData{i}(43:48));%dep
                uOutput(i,7) = str2num(mData{i}(44:51));%dep
            end

            %             uOutput(i,8) = str2num(mData{i}(12:13));    %hr
            %             uOutput(i,9) = str2num(mData{i}(15:16));    %min
            uOutput(i,8) = str2num(mData{i}(14:15));    %hr
            uOutput(i,9) = str2num(mData{i}(16:17));    %min
            %uOutput(i,:);
        catch
            msg.dbfprintf('Import: Problem in line %d of %s. Line ignored.\n',i, sFilename);
            uOutput(i,:)=nan;
        end
    end
    l = isnan(uOutput(:,1));
    uOutput(l,:) = [];
end

author: CelsoReyes


text:
% ------------------------------------------------------
% SwarmOps - Heuristic optimization for Matlab
% Copyright (C) 2003-2010 Magnus Erik Hvass Pedersen.
% Please see the file license.txt for license details.
% SwarmOps on the internet: http://www.Hvass-Labs.org/
% ------------------------------------------------------

% Many Optimizing Liaisons (MOL) optimization method is a
% simplification of PSO originally by Eberhart et al. (1, 2).
% MOL does not have any attraction to the particle's own best
% known position. It is similar to the "Social Only" PSO
% suggested by Kennedy (3) and was studied more thoroguhly
% by Pedersen et al. (4) who found it to sometimes outperform
% PSO and have behavioural parameters that were easier to tune.
% Literature references:
% (1) J. Kennedy and R. Eberhart. Particle swarm optimization.
%     In Proceedings of IEEE International Conference on Neural
%     Networks, volume IV, pages 1942-1948, Perth, Australia, 1995
% (2) Y. Shi and R.C. Eberhart. A modified particle swarm optimizer.
%     In Proceedings of the IEEE International Conference on
%     Evolutionary Computation, pages 69-73, Anchorage, AK, USA, 1998.
% (3) J. Kennedy. The particle swarm: social adaptation of knowledge,
%     In: Proceedings of the IEEE International Conference on
%     Evolutionary Computation, Indianapolis, USA, 1997.
% (4) M.E.H. Pedersen and A.J. Chipperfield. Simplifying particle swarm
%     optimization. Applied Soft Computing, volume 10, pages 618-628, 2010.
% Parameters:
%     problem; name or handle of optimization problem, e.g. @myproblem.
%     data; data-struct, see e.g. the file myproblemdata.m
%     parameters; behavioural parameters for optimizer,
%                 see file molparameters.m
% Returns:
%     bestX; best found position in the search-space.
%     bestFitness; fitness of bestX.
%     evaluations; number of fitness evaluations performed.
function [bestX, bestFitness, evaluations] = mol(problem, data, parameters)

    % Copy data contents to local variables for convenience.
    n = data.Dim;
    acceptableFitness = data.AcceptableFitness;
    maxEvaluations = data.MaxEvaluations;
    lowerBound = data.LowerBound;
    upperBound = data.UpperBound;

    % Behavioural parameters for this optimizer.
    s = parameters(1);        % Swarm-size
    omega = parameters(2);    % Inertia weight.
    phiG = parameters(3);     % Swarm's best weight.

    % Initialize the velocity boundaries.
    range = upperBound-lowerBound;
    lowerVelocity = -range;
    upperVelocity = range;

    % Initialize swarm.
    x = initpopulation(s, n, data.LowerInit, data.UpperInit);    % Particle positions.
    v = initpopulation(s, n, lowerVelocity, upperVelocity);      % Velocities.

    % Determine swarm's best-known position and its fitness.
    bestFitness = Inf;
    for i=1:s
        % Compute fitness for i'th agent.
        fitness = feval(problem, x(i,:), data);

        % If improvement, update best-known position and fitness.
        if fitness < bestFitness
            bestFitness = fitness;
            bestX = x(i,:);
        end
    end

    % Perform optimization iterations until acceptable fitness
    % is achieved or the maximum number of fitness evaluations
    % has been performed.
    evaluations = s; % Fitness evaluations above count as iterations.
    while (evaluations < maxEvaluations) && (bestFitness > acceptableFitness)

        % Pick index for a random particle from the swarm.
        i = randindex(s);

        % Pick random weight.
        rG = rand(1, 1);

        % Update velocity for i'th particle.
        v(i,:) = omega * v(i,:) + rG * phiG * (bestX - x(i,:));

        % Bound velocity.
        v(i,:) = bound(v(i,:), lowerVelocity, upperVelocity);

        % Update position for i'th particle.
        x(i,:) = x(i,:) + v(i,:);

        % Bound position to search-space.
        x(i,:) = bound(x(i,:), lowerBound, upperBound);

        % Compute fitness.
        fitness = feval(problem, x(i,:), data);

        % Update swarm's best-known position.
        if fitness < bestFitness
            % Update swarm's best-known fitness.
            bestFitness = fitness;

            % Update swarm's best-known position.
            % This must be copied because the particles
            % will continue to move in the search-space.
            bestX = x(i,:);
        end

        % Increment counter.
        evaluations = evaluations + 1;
    end
end

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


text:
classdef ClosestVigdergauzSuperEllipseComputer < handle

    properties (Access = public)
        xopt
        error
    end

    properties (Access = private)
        txi
        rho
        q
        problem
        frames
        optimalExponent
    end

    methods (Access = public)

        function  obj = ClosestVigdergauzSuperEllipseComputer()
            obj.init();
            obj.solveProblem();
            obj.writeOptimizationVideo();
            obj.printVigdergauzSuperEllipse();
        end

    end

    methods (Access = private)

        function init(obj)
            obj.txi = pi/8;
            obj.rho = 0.8;
        end

        function solveProblem(obj)
            s.rho = obj.rho;
            s.txi = obj.txi;
            s.savingFrames = true;
            obj.optimalExponent = OptimalExponentClosestToVigergauz(s);
            obj.optimalExponent.compute();
            obj.xopt   = obj.optimalExponent.qOpt;
            obj.q      = obj.optimalExponent.qOpt;
            obj.error  = obj.optimalExponent.error;
            obj.frames = obj.optimalExponent.frames;
        end

        function writeOptimizationVideo(obj)
            v = VideoWriter('SuperEllipseToVigdergauz.avi');
            open(v);
            nf = 6;
            for iframe = 1:nf*(numel(obj.frames)-2)
                i = floor(iframe/nf)+1;
                frame = obj.frames{i};
                writeVideo(v,frame);
            end
            close(v);
        end

        function printVigdergauzSuperEllipse(obj)
            s.mx = obj.optimalExponent.mx;
            s.my = obj.optimalExponent.my;
            s.txi = obj.txi;
            s.rho = obj.rho;
            s.q   = obj.q;
            stressPrinter = VigdergauzSuperEllipseStressPrinter(s);
            stressPrinter.print();
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [Angles, Offsets, w1cw, w1rms, w1rp, Tau] = SPGR_prepare( Prot )

Angles  = Prot.Angles;
Offsets = Prot.Offsets;
w1cw  = zeros(length(Angles),1);
w1rms = w1cw;
w1rp  = w1cw;
Tau   = w1cw;

for ii = 1:length(Angles)
    Pulse = GetPulse(Angles(ii), Offsets(ii), Prot.Tm, Prot.MTpulse.shape, Prot.MTpulse.opt);
    w1cw(ii) = compute_w1cw(Prot.TR, Pulse);
    w1rms(ii) = compute_w1rms(Pulse);
    [w1rp(ii), Tau(ii)] = compute_w1rp(Pulse);
end

end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function arrivalobj = singlestation_specrat(dbpath, expr, f1, f2, pretrigger, posttrigger, max_arrivals)
%HANKELQ Compute spectral ratio in two frequency bands
%   HANKELQ(dbpath, expr, freq_high, freq_low, pretrigger_seconds, posttrigger_seconds, max_arrivals )
%   Loads arrivals from an arrival table (after subsetting with expr)
%   retrieves waveform data corresponding to those arrivals, cleans and
%   plots the waveform data. The spectral ratio in bands around freq_high
%   and freq_low is computed. When done for multiple earthquakes, a plot of this the
%   natural log of this ratio versus travel time has a slope from which Q
%   can be measured.
%   pretrigger_seconds is the number of seconds before the arrival time to
%   get waveform data for.
%   posttrigger_seconds is the number of seconds after the arrival time to
%   get waveform data for.
%   max_arrivals is the maximum number of arrivals to process.
%
%     Based on the method of Arthur Hankel, "The Effects of Attenuation and
%     Site Response on the Spectra of Microearthquakes in the Northeastern
%     Caribbean", BSSA, 72, 4, 1379-1402.
%
%   Example:
%
%
%   History:
%     April 2014: Glenn Thompson
%       Original version: load arrivals, load waveforms and plot waveforms
%     April-November 2014: Heather McFarlin
%       Modified to also plot amplitude spectra
%     November 20, 2014: Glenn Thompson
%       Completely overhauled & modularized.
%       Modified method to compute amplitude spectra. Now computes Q too. Now
%       also generic - works with input variables so it can be used on different
%       databases, for example.
%     November 2017: Glenn Thompson
%       Fixed to work with updated GISMO classes
%                    Heather McFarlin
%       Added units to Y-axis on plots
%       Added amplitude spectrum of noise before waveform onto amplitude spectrum plot
%       Added noise part of waveform to top subplot
%     NEED TO FIX TO CALCULATE Q based on seaz- Q is azimuthally dependent
%     at Uturuncu stations!

    if ~admin.antelope_exists
        warning('Antelope not installed on this computer')
        return
    end

    if ~exist('max_arrivals','var')
        max_arrivals=Inf;
    end

    taper_seconds=pretrigger+posttrigger;

    %arrivals = antelope.dbgetarrivals(dbpath, expr);
    arrivalobj = Arrival.retrieve('antelope', dbpath, 'subset_expr', expr);
    %arrivalobj = arrivalobj.subset(expr)
    %arrivalobj = arrivalobj.subset(1:max_arrivals);  %when done, comment this line out, it's just for testing
    arrivalobj = arrivalobj.addwaveforms(datasource('antelope', dbpath), pretrigger+taper_seconds, posttrigger+taper_seconds);

%     w = antelope.arrivals2waveforms(dbpath, arrivals, pretrigger, posttrigger, taper_seconds, max_arrivals);
    %w = waveform_clean(w);
    close all

    %[y, t]=plot_arrival_waveforms(arrivals, w, pretrigger, posttrigger, taper_seconds, max_arrivals, f1, f2);
    [y, t]=plot_arrival_waveforms(arrivalobj, pretrigger, posttrigger, taper_seconds, max_arrivals, f1, f2);
    % This is the figure we use to derive q from spectral ratios for
    % earthquakes of different distances (travel times)
    % Q comes from the slope
    y
    t
    figure;
    plot(t, y, 'o');
    ylabel('ln(A_1/A_2)')
    xlabel('travel time (s)')
    p = polyfit(t,y,1);
    xlim=get(gca,'XLim');
    yline = polyval(p, xlim);
    hold on;
    plot(xlim, yline)

    slope = (yline(2) - yline(1)) / (xlim(2) - xlim(1));
    q = - pi * (f1 - f2) / slope;
    title(sprintf('Q = %.0f', q)); %% Can we add something that denotes whether the Q value is for Qp or Qs?
end

%function [y, t] = plot_arrival_waveforms(arrivals, w, pretrigger, posttrigger, taper_seconds, max_arrivals, f1, f2)
function [y, t] = plot_arrival_waveforms(arrivalobj, pretrigger, posttrigger, taper_seconds, max_arrivals, f1, f2)

    FMIN = 1.0;

    %% open an output file
    fid=fopen([mfilename,'.txt'],'w');

    taper_fraction = (taper_seconds * 2) / (taper_seconds * 2 + pretrigger + posttrigger);

    % get travel time for p-wave
    p_time = arrivalobj.time - arrivalobj.otime
    anum = arrivalobj.time;

    w = arrivalobj.waveforms

    for i=1:min([max_arrivals numel(w)])
        signal = get(w(i),'data');
        N = length(signal);
        if N>0

            % taper, filter to remove long period noise, and cut out tapered signal to leave only
            % untapered signal
            wf = w(i);

            taperwin=tukeywin(N, taper_fraction);
            signal=signal.*taperwin;
            wf = set(wf, 'data', signal);
            fmax = get(wf, 'freq') * 0.4;
            fobj = filterobject('b', [FMIN fmax], 2);
            wf = filtfilt(fobj, wf);
            [snum, enum]=gettimerange(wf);

            %get noise before event
            wf_noise = extract(wf, 'time', snum, anum(i)-pretrigger/86400)

            %wf=subtime(wf, snum+taper_seconds/86400, enum-taper_seconds/86400);
            wf=extract(wf, 'time', snum+taper_seconds/86400, enum-taper_seconds/86400)

            % integrate
            dwf = integrate(wf);

            % integrate noise

            dwf_noise = integrate(wf_noise);

            % Plot waveform
            hf(i)=figure;
            ax(i,1)=subplot(3,1,1);
            plot(w(i), 'xunit', 'date', 'color', 'r', 'axeshandle', ax(i,1)); %detrended and cleaned waveform -artifacts at beginning of signal are likely due to cleaning and detrending filters
            datetick('x','HH:MM:SS', 'keeplimits')
            hold on;
            plot(wf, 'xunit', 'date', 'color', 'g', 'axeshandle', ax(i,1)); %filtered signal
            datetick('x','HH:MM:SS','keeplimits')
            hold on
            plot(wf_noise, 'xunit', 'date', 'color', 'k', 'axeshandle', ax(i,1)) % noise
            datetick('x','HH:MM:SS', 'keeplimits')
            xlabel(sprintf('Time with arrival at %s',datestr(anum(i), 'yyyy-mm-dd HH:MM:SS.FFF')));
            ylabel('Velocity (nm/s)');
            title(sprintf('%s',arrivalobj.channelinfo{i}));
            % plot arrival time as grey dotted line
            ylim=get(gca,'YLim');
            hold on
            plot([anum(i) anum(i)],ylim,'Color',[0.5 0.5 0.5], 'linestyle', '--')
            hold off

            ax(i,2)=subplot(3,1,2);
            plot(dwf, 'xunit', 'date', 'color', 'b', 'axeshandle', ax(i,2)); %filtered & integrated signal
            datetick('x','HH:MM:SS:FFF', 'keeplimits', 'keepticks')
            xlabel(sprintf('Time with arrival at %s',datestr(anum(i), 'yyyy-mm-dd HH:MM:SS.FFF')));
            ylabel('Displacement (nm)');
            title(sprintf('%s',arrivalobj.channelinfo{i}));
            % plot arrival time as grey dotted line
            ylim=get(gca,'YLim');
            hold on
            plot([anum(i) anum(i)],ylim,'Color',[0.5 0.5 0.5], 'linestyle', '--')
            hold off

            % compute and plot amplitude spectrum
            s = amplitude_spectrum(dwf);
            A = s.amp;
            f = s.f;
            phi = s.phi;
%             [A, phi, f] = amplitude_spectrum(dwf);
            ax(i,3)=subplot(3,1,3);
            A=smooth(A);
            semilogy(f,A);
            size(f)
            %loglog(f, A)
            xlabel('Frequency (Hz)')
            ylabel('Amplitude (nm)')
            hold on
            % add noise spectrum
            s_noise = amplitude_spectrum(dwf_noise);
            A_noise = s_noise.amp;
            f_noise = s_noise.f;
            phi_noise = s_noise.phi;
            A_noise = smooth(A_noise);
            semilogy(f_noise, A_noise, 'color', 'k');
            size(f_noise)
            hold on
            % evaluate the spectral ratio
            A1=mean(A(find(f>=f1*0.9 & f<=f1*1.1))); %  for 20 Hz, this is 18-22 Hz
            A2=mean(A(find(f>=f2*0.8 & f<=f2*1.2 ))); % for 5 Hz this is 4-6 Hz
            patch([f1*0.9 f1*1.1 f1*1.1 f1*0.9],[0 0 A1 A1],[0.9 0.9 0.9]);
            patch([f2*0.8 f2*1.2 f2*1.2 f2*0.8],[0 0 A2 A2],[0.9 0.9 0.9]);
            %evaluate spectral ratio of  noise
            A1_noise = mean(A_noise(find(f_noise>=f1*0.9 & f_noise<=f1*1.1))); % for 20 Hz noise
            A2_noise = mean(A_noise(find(f_noise>= f2*0.8 & f_noise<=f2*1.2))); % for 5 Hz noise
            snr1 = A1/A1_noise;
            snr2 = A2/A2_noise;
            disp(sprintf('A1 = %6.3f, A2 = %6.3f, A1/A2 = %5.2f', A1, A2, A1/A2));
            hold on
            disp(sprintf('SNR1 = %6.3f, SNR2 = %6.3f', snr1, snr2));
            hold off

            % add title
            %outstr=sprintf('phase=%s orid=%d seaz=%6.2f delta=%5.3f depth=%5.1f\nA1=%5.3f A2=%5.3f A1/A2=%5.2f\n',arrivals.phase{i}, arrivals.orid(i), arrivals.seaz(i), arrivals.delta(i), arrivals.depth(i), A1, A2, A1/A2);
            outstr=sprintf('phase=%s orid=%d seaz=%6.2f delta=%5.3f depth=%5.1f\nA1=%5.3f A2=%5.3f A1/A2=%5.2f\n SNR1 = %6.3f SNR2 = %6.3f',arrivalobj.iphase{i}, arrivalobj.orid(i), arrivalobj.seaz(i), arrivalobj.delta(i), arrivalobj.depth(i), A1, A2, A1/A2, snr1, snr2);
            title(outstr)

            % write out to file & close plot
            filename=sprintf('%s-%d',mfilename,i);
            print('-dpng', figure(i),filename)
            if numel(w)>100  % I changed this from 10 to 100 because when I am passing more than 10 arrivals, my figures end up blank
                close
            end

            % compute y=ln A1/A2 & t for formula 2 in Frankel (1982)
            y(i) = log(A1/A2);
            [times, phasenames] = arrtimes(arrivalobj.delta(i), arrivalobj.depth(i));
            phase_index = 1;
            found = false;
            while phase_index <= length(times) & ~found
                thisphase = lower(phasenames{phase_index});
                thisphase = thisphase(1);
                %if strcmp(lower(arrivals.phase{i}), thisphase)
                if strcmp(lower(arrivalobj.iphase{i}), thisphase)
                    t(i)=times(phase_index);
                    found=true;
                end
                phase_index = phase_index + 1;
            end

            %% write out to file
            fprintf(fid,'%14.2f %8.2f %6.3f %5.1f %4.1f %4.1f %e %e\n', p_time(i), t(i), arrivalobj.delta(i), arrivalobj.depth(i), f1, f2, A1, A2);

        end

    end

    %% close output file
    fclose(fid);
end

author: geoscience-community-codes

repo: GISMO


text:
%% This method normalizes the features to fall in the range [-1, 1]

function testX = testFeatNormalize(testX, minimums, ranges)

	testX = (testX - repmat(minimums, size(testX, 1), 1)) ./ repmat(ranges, size(testX, 1), 1);

end

author: vztu

repo: VIDEVAL

save_path: github-repos/MATLAB/vztu-VIDEVAL


text:
function epo= proc_selectEpochs(epo, idx, varargin)
%PROC_SELECTEPOCHS - Select specified epochs within epoch structure
%
%Synopsis:
%  EPO= proc_selectEpochs(EPO, IDX, <OPT>)
%  EPO= proc_selectEpochs(EPO, 'not', IDX, <OPT>)
%
%Arguments:
%  EPO: STRUCT - Epoch structure
% selects the epochs with specified indices. void classes are removed.
%
% the structure 'epo' may contain a field 'indexedByEpochs' being a
% cell array of field names of epo. in this case subarrays of those
% fields are selected. here it is assumed that the last dimension
% is indexed by epochs.
%
% IN  epo  - structure of epoched data
%     idx  - indices of epochs that are to be selected [or are to be
%            excluded if the 'not' form is used].
%            if this argument is omitted, all non-rejected epochs are
%            selected, i.e., epochs with any(epo.y).
%
% OUT epo  - updated data structure

% Benjamin Blankertz
props= {'RemoveVoidClasses',   1,   'BOOL'};

if nargin==0,
  epo= props; return
end

misc_checkType(epo, 'STRUCT(x clab)');
epo= misc_history(epo);

if isequal(idx, 'not'),
  idx= setdiff(1:size(epo.y,2), varargin{1},'legacy');
  varargin= varargin(2:end);
end

opt= opt_proplistToStruct(varargin{:});
opt= opt_setDefaults(opt, props);

if ~exist('idx','var'),
  %% select accepted epochs
  idx= find(any(epo.y==1,1));
end

subidx= cat(2, repmat({':'}, 1, ndims(epo.x)-1), {idx});
epo.x= epo.x(subidx{:});

if isfield(epo, 'y'),
  epo.y= epo.y(:,idx);

  nonvoidClasses= find(any(epo.y==1,2));
  if length(nonvoidClasses)<size(epo.y,1) && opt.RemoveVoidClasses
    msg= sprintf('void classes removed, %d classes remaining', ...
                 length(nonvoidClasses));
    warning(msg, 'selection', mfilename);
    epo.y= epo.y(nonvoidClasses,:);
    if isfield(epo, 'className'),
      epo.className= {epo.className{nonvoidClasses}};
    end
  end
end

if isfield(epo, 'event'),
  for Fld= fieldnames(epo.event)',
    fld= Fld{1};
    tmp= getfield(epo.event, fld);
    % the first dimension must be indexed by epochs
    subidx= repmat({':'}, 1, ndims(tmp));
    subidx{1}= idx;
    epo.event= setfield(epo.event, fld, tmp(subidx{:}));
  end
end

author: bbci


text:
function [] = PUMA_InverseKinematics(Trans, theta, th, dstPoint)
%% Variable Initialization

% Variables to store error related data
Kp = 0.0045;
errorThreshold = 0.1;

% Matrix containing Joint angle Information and Desired Point Location
t1 = th(1); t2 = th(2); t3 = th(3);
Theta = [double(subs(theta(1))); double(subs(theta(2))); double(subs(theta(3)))];

% Make a column matrix representing the Position of the End-Effector
endPoint = [Trans(1,4,6); Trans(2,4,6); Trans(3,4,6)];
% Find the current numeric end-effector position
currPoint = double(subs(endPoint));
% Find the error vector
errorVec = (dstPoint - currPoint);
% Create a Jacobian Matrix using the given Theta values
Jacobian = PUMA_Jacob(Theta(1,1), Theta(2,1), Theta(3,1));

for i=1:3
    [X(i), Y(i), Z(i)] = PUMA_Plot(double(subs(Trans(:,:,i))));
end
for i=4:6
    X(i) = double(subs(Trans(1,4,i)));
    Y(i) = double(subs(Trans(2,4,i)));
    Z(i) = double(subs(Trans(3,4,i)));
end

PUMA_Draw(X, Y, Z);

%% Inverse Kinematics Initialization

% Apply the control equation to the error vector
deltaTheta = transpose(Jacobian) * errorVec * Kp;
% Apply the control changes to the Joint Angles
Theta(1,1) = Theta(1,1) + deltaTheta(1,1);
Theta(2,1) = Theta(2,1) + deltaTheta(2,1);
Theta(3,1) = Theta(3,1) + deltaTheta(3,1);
t1 = Theta(1,1); t2 = Theta(2,1); t3 = Theta(3,1);

% Find the euclidean distance between the desired and current position of
% the manipulator.
dist = sqrt((dstPoint(1,1) - currPoint(1,1))^2 + ...
    (dstPoint(2,1)-currPoint(2,1))^2 +(dstPoint(3,1)-currPoint(3,1))^2);

%% Iterative Inverse Kinematics for correcting the manipulator position.

% Correct the manipulator position to reduce the distance.
while (dist > errorThreshold)
    Jacobian = PUMA_Jacob(Theta(1,1), Theta(2,1), Theta(3,1));
    currPoint = double(subs(endPoint));
    errorVec = (dstPoint - currPoint);
    deltaTheta = transpose(Jacobian) * errorVec * Kp;
    Theta(1,1) = Theta(1,1) + deltaTheta(1,1);
    Theta(2,1) = Theta(2,1) + deltaTheta(2,1);
    Theta(3,1) = Theta(3,1) + deltaTheta(3,1);
    t1 = Theta(1,1); t2 = Theta(2,1); t3 = Theta(3,1);
    dist = sqrt((dstPoint(1,1) - currPoint(1,1))^2 + ...
        (dstPoint(2,1)-currPoint(2,1))^2 +(dstPoint(3,1)-currPoint(3,1))^2);

    for i=1:3
        [X(i), Y(i), Z(i)] = PUMA_Plot(double(subs(Trans(:,:,i))));
    end
    for i=4:6
        X(i) = double(subs(Trans(1,4,i)));
        Y(i) = double(subs(Trans(2,4,i))); Z(i) = double(subs(Trans(3,4,i)));
    end
    PUMA_Draw(X, Y, Z);

end

cla
for i=1:3
    [X(i), Y(i), Z(i)] = PUMA_Plot(double(subs(Trans(:,:,i))));
end
for i=4:6
    X(i) = double(subs(Trans(1,4,i)));
    Y(i) = double(subs(Trans(2,4,i)));
    Z(i) = double(subs(Trans(3,4,i)));
end
PUMA_Draw(X, Y, Z);
end

author: YashBansod

repo: Robotics-Planning-Dynamics-and-Control


function quer(str);
% function quer(str);
% Query user
% if str=='c' Continue or not
% if str=='g' close Graphs or not
str=upper(str);
switch str
case 'C'

    button=questdlg('Continue? Y/N',...
        'Continue Program',...
        'Yes','No','Yes');
    switch button
    case 'No'
        error('Program execution stopped');
    case 'Yes'
        fprintf('%-30%s \n','  ');
    end
case 'G'
    button=questdlg('Close All Graphs Y/N?','Graphs','No','Yes','Yes');
    switch button
    case 'No'
        fprintf('%-30%s \n','  ');
    case 'Yes'
        close all;
    end
end


text:
function F = tan(F, varargin)
%TAN   Tangent of a CHEBFUN.
%   TAN(F) computes the tangent of the CHEBFUN F.
%
%   TAN(F, PREF) does the same but uses the CHEBFUNPREF object PREF when
%   computing the composition.
%
% See also ATAN, TAND.

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

% [TODO]:  Restore or change this once we have decided the proper behavior or
% isfinite() and defined that function.
% if ( ~isfinite(f) )
%     error('CHEBFUN:CHEBFUN:tan:inf',...
%         'SIN is not defined for functions which diverge to infinity');
% end

% Call the compose method:
F = compose(F, @tan, varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function loss_group_info=gen_network_loss_group(train_opts, net_config, group_input_info)

feat_output_dim=group_input_info.var_dims;
assert(length(feat_output_dim)==1);

loss_group_info=My_net_util.gen_group_info_basic(net_config, 'loss_group');

stageout_info=[];
stageout_info.stage_idx=1;
stageout_info.output_dim=feat_output_dim;
stageout_info.stageout_num=1;
stageout_info.gen_prediction_info=true;

[loss_group_info.net_info, dag_net_lossgroup]=do_gen_net_info(train_opts, stageout_info);

if ~isempty(dag_net_lossgroup)
    loss_group_info.dag_net=dag_net_lossgroup;
    loss_group_info.use_dagnn=true;
end

loss_group_info.forward_evaluate_fn=@cnn_loss_group_evaluate;

loss_group_info.prediction_layer_idxes=loss_group_info.net_info.ref.prediction_layer_idxes;
net_config.ref.group_infos{loss_group_info.group_idx,1}=loss_group_info;

net_config.ref.root_group_info.child_group_idxes(end+1,1)=loss_group_info.group_idx;

end

function [net_info, dag_net]=do_gen_net_info(train_opts, stageout_info)

    loss_config=train_opts.loss_config;
    class_num=loss_config.class_num;
    loss_group_input_dim=stageout_info.output_dim;

    layers=cell(0);

    dag_net=[];
    if loss_config.lossgroup_conv_num>0
        [dag_net, output_dim_crf_conv1, one_layers]=gen_dag_net_lossgroup(loss_config, loss_group_input_dim);
        layers=cat(2, layers, one_layers);
    else
        output_dim_crf_conv1=loss_group_input_dim;
    end

    layers{end+1} = struct('type', 'dropout', 'rate', 0.5) ;

    output_dim_crf_conv2=class_num;
    input_dim_crf_conv2=output_dim_crf_conv1;
    layers{end+1} = struct('type', 'conv', ...
                   'filters', 0.01 * randn(3,3,input_dim_crf_conv2, output_dim_crf_conv2,'single'), ...
                   'biases', zeros(1, output_dim_crf_conv2, 'single'), ...
                   'stride', 1, ...
                   'pad', 0 );

    one_layer=[];
    one_layer.type='my_custom';
    one_layer.custom_type='dense_softmaxloss';

    one_layer.forward_fn=@cnn_layer_dense_softmaxloss_forward;
    one_layer.backward_fn=@cnn_layer_dense_softmaxloss_backward;

    one_layer.stageout_info=stageout_info;

    one_layer.layer_update_fn=[];
    layers{end+1}=one_layer;

layers=gen_padding_keep_size(layers);

net_info=My_net_util.gen_net_info_basic();

net_info.ref.name='final_loss';
net_info.ref.layers=layers;
net_info.ref.prediction_layer_idxes=length(layers);

net_info.ref.lr_multiplier=loss_config.lr_multiplier;

end

function [dag_net, output_dim, layers]=gen_dag_net_lossgroup(loss_config, loss_group_input_dim)

conv_num=loss_config.lossgroup_conv_num;
one_output_dim=loss_config.lossgroup_conv_filter_num;

dag_net=dagnn.DagNN();

start_var_name='data_input';

layer_gen_info=[];
layer_gen_info.one_output_dim=loss_group_input_dim;
layer_gen_info.one_outputs={start_var_name};

if loss_group_input_dim~=one_output_dim
	error('should not come here!');
end

layer_name_prefix=sprintf('lossblock');
layer_gen_info=My_net_util.add_res_conv_block(dag_net, layer_gen_info, one_output_dim, conv_num, layer_name_prefix);

output_dim=layer_gen_info.one_output_dim;

output_var_name=layer_gen_info.one_outputs{1};
output_var_idx=dag_net.getVarIndex(output_var_name);
dag_net.vars(output_var_idx).fanout=0;

% debug:
% dag_net.print('Format', 'dot');

My_net_util.fix_padding_resnet(dag_net);

layers=cell(0);

one_layer=[];
one_layer.type='my_custom';
one_layer.custom_type='dagnn_wrapper';
one_layer.forward_fn=@cnn_layer_dagnn_wrapper_forward;
one_layer.backward_fn=@cnn_layer_dagnn_wrapper_backward;
one_layer.layer_update_fn=[];

one_layer.input_var_names=dag_net.getInputs();
one_layer.output_var_names=dag_net.getOutputs();

assert(length(one_layer.input_var_names)==1);
assert(length(one_layer.output_var_names)==1);

one_layer.use_single_output=true;

layers{end+1}=one_layer;

end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function [stateLog, errorStr, errorEventNum] = ma_produceStateLogFromScript(script, maData, celBodyData, allowInterrupt, varargin)
%ma_produceStateLogFromScript Summary of this function goes here
%   Detailed explanation goes here
    if(isempty(varargin))
        stateLog = [];
        eventNum = 1;
    else
        partialExec = varargin{1};
        stateLog = partialExec{1};
        eventNum = partialExec{2};
    end
    initialEventNum = eventNum;

    if(initialEventNum > length(script))
        initialEventNum = length(script);
    end

%     profile on -history
    errorStr = '';
    errorEventNum = -1;
    for(i = initialEventNum:length(script)) %#ok<*NO4LP>
        event = script{i};

        if(allowInterrupt && mod(i,2))
            drawnow;
        end

        try
            eventLog = ma_executeEvent(event, getLastStateLogEntry(stateLog), eventNum, maData, celBodyData);
            stateLog = [stateLog ; eventLog]; %#ok<*AGROW>
            eventNum = eventNum + 1;
        catch ME
            errorStr = ME.message;
            errorEventNum = eventNum;

            disp('################## Propagation Error ###################');
            disp(['Event Num: ', num2str(eventNum)]);
            disp(['Error Msg: ', ME.message]);
            disp('#########################################################');
            try
                disp(ME.cause{1}.message);
                disp('#########################################################');
            catch
            end
            for(i=1:length(ME.stack)) %#ok<FXSET>
                disp(['Index: ', num2str(i)]);
                disp(['File: ',ME.stack(i).file]);
                disp(['Name: ',ME.stack(i).name]);
                disp(['Line: ',num2str(ME.stack(i).line)]);
                disp('####################');
            end
            disp('#########################################################');
            break;
        end
    end
%     profile viewer

    stateLog(all(stateLog==0,2),:)=[];
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef DehomogenizingRadialExample < handle

    properties (Access = private)
        backgroundMesh
        theta
        superEllipse
        cellLevelSetParams
    end

    properties (Access = private)
        nx1
        nx2
        nCells
    end

    methods (Access = public)

        function obj = DehomogenizingRadialExample()
            obj.init();
            for i = 10:50
                obj.nCells = i;
                obj.createBackgroundMesh();
                obj.createOrientation();
                obj.createSuperEllipseParams();
                obj.createLevelSetCellParams();
                obj.dehomogenize();
                exportgraphics(gcf,'testAnimated.gif','Append',true);
            end
        end

    end

    methods (Access = private)

        function init(obj)
            obj.nx1    = 40*2;%180
            obj.nx2    = 40;%180
            obj.nCells = 32;%32
        end

        function createBackgroundMesh(obj)

            x1 = linspace(-1,1,obj.nx1);
            x2 = linspace(0,1,obj.nx2);

            x1T = repmat(x1,obj.nx2,1);
            x2T = repmat(x2',1,obj.nx1);

%             xy = obj.coord;
%             x1 = xy(:,1);
%             x2 = xy(:,2);
%             s.connec   = delaunay(x1,x2);
%             s.coord    = obj.coord;
%             obj.backgroundMesh = Mesh(s);

%             x1min = min(x1);
%             x1max = max(x1);
%             x2min = min(x2);
%             x2max = max(x2);
%             [coordinates, nodes,nel,nnode] = MeshRectangular(x1max-x1min,x2max-x2min,obj.nx1,obj.nx2);
%             s.coord(:,1) = coordinates(:,1)+x1min;
%             s.coord(:,2) = coordinates(:,2)+x2min;
%             s.connec = nodes;
%             obj.backgroundMesh = Mesh(s);
%             obj.backgroundMesh.plot();

             [xv,yv] = meshgrid(x1,x2);
             [F,V] = mesh2tri(xv,yv,zeros(size(xv)),'x');
             s.coord  = V(:,1:2);
             s.connec = F;
             obj.backgroundMesh = Mesh(s);
        %     obj.backgroundMesh.plot()
%             obj.coord = s.coord;

        end

        function createOrientation(obj)
            x2 = obj.backgroundMesh.coord(:,2);
            x1 = obj.backgroundMesh.coord(:,1);
            obj.theta = atan((x2+0.1)./x1);
            isLeft = x1 < 0;
            %obj.theta(isLeft) = obj.theta(isLeft) + 180;
        end

        function createSuperEllipseParams(obj)
           s.coord = obj.backgroundMesh.coord;
           s.mMin  = 0.4;
           s.mMax  = 0.99;
           s.qMin  = 32;
           s.qMax  = 32;
           sE = SuperEllipseDistributionExample(s);
           sE.computeParameters();
           obj.superEllipse = sE;
        end

        function createLevelSetCellParams(obj)
           s.type   = 'rectangleInclusion';%'smoothRectangle';
           s.widthH = 0.87*ones(size(obj.superEllipse.m1));
           s.widthV = 0.87*ones(size(obj.superEllipse.m2));
           s.pnorm  = obj.superEllipse.q;
           s.ndim   = 2;
           obj.cellLevelSetParams = s;
        end

        function dehomogenize(obj)
            s.backgroundMesh     = obj.backgroundMesh;
            s.nCells             = obj.nCells;
            s.theta              = obj.theta;
            s.cellLevelSetParams = obj.cellLevelSetParams;
            s.mesh               = obj.backgroundMesh;
            d = Dehomogenizer(s);
            d.compute();
            d.plot();
        end

    end

end


text:
function close(freader)
%
%   O = close(freader) - closes the binary file
%
fclose(freader.fd);
freader.fd = 0;

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function  eta = conftim1(a, om, omr, omt, h)
% calculate the conformal time today for a universe with a cosmological constant
% oml = ratio of dark energy and matter = ohm_l/ohm_m
%
% D Vangheluwe 2 oct 2004
% modified for curvature of space with extra variable : omt (dd 31 mrt 2005)

global  GL_cmb_c  GL_cmb_h0;

% velocity of light in m/s
%c = 2.998e8;
c = GL_cmb_c;
% the hubble constant at present in h Mpc^-1, see my notes p71
%h0 = 1e5/c;
h0 = GL_cmb_h0;

eta = 1 ./ (h0 * h * sqrt(omr + om * a + (omt - omr - om) * a.^4 + (1 - omt) * a.^2));

author: Sable


text:
%% Summarize Text Using Transformers
% This example shows how to summarize a piece of text using GPT-2.
%
% Transformer networks such as GPT-2 can be used to summarize a piece of
% text. The trained GPT-2 transformer can generate text given an initial
% sequence of words as input. The model was trained on comments left on
% various web pages and internet forums.
%
% Because lots of these comments themselves contain a summary indicated by
% the statement "TL;DR" (Too long, didn't read), you can use the
% transformer model to generate a summary by appending "TL;DR" to the input
% text. The |generateSummary| function takes the input text, automatically
% appends the string |"TL;DR"| and generates the summary.

%% Load Transformer Model
% Load the GPT-2 transformer model using the |gpt2| function.

mdl = gpt2;

%% Load Data
% Extract the help text for the |eigs| function.

inputText = help('eigs')

%% Generate Summary
% Summarize the text using the |generateSummary| function.

rng('default')
summary = generateSummary(mdl,inputText)

author: matlab-deep-learning

repo: transformer-models


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

%*****************************************************************************80
%
%% RULE03 returns the rule of degree 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 July 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, integer N, the number of nodes.
%
%    Output, real X(2,N), the coordinates of the nodes.
%
%    Output, real W(N), the weights.
%
  xs = [ ...
    -.5773502691896256,0.5773502691896260, ...
    0.5773502691896256,-.5773502691896260 ];
  ys = [ ...
    -.5773502691896260,-.5773502691896256, ...
    0.5773502691896260,0.5773502691896256 ];
  ws = [ ...
    0.7071067811865476,0.7071067811865476, ...
    0.7071067811865476,0.7071067811865476 ];

  x(1,1:n) = xs(1:n);
  x(2,1:n) = ys(1:n);
  w(1:n) = ws(1:n);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function pass = test_CLA( pref )
% Try out some basic continuous linear algebra.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 200*pref.cheb2Prefs.chebfun2eps;

seedRNG(0)
gam = 10;
% Add some Gaussians:
f = @(x, y) 0;
for n = 1:20
    x0 = 2*rand-1;
    y0 = 2*rand-1;
    df = @(x,y) exp(-gam*((x-x0).^2 + (y-y0).^2));
    f = @(x, y) f(x, y) + df(x, y);
end

% Make a CHEBFUN2:
F = chebfun2(f);

% Check QR of a chebfun2

% Test QR of the columns:
[Q, R] = qr(F.cols);
pass(1) = ( norm(Q*R - F.cols) < tol );

% Test QR of the rows:
[Q, R] = qr(F.rows);
pass(2) = ( norm(Q*R - F.rows) < tol );

% Check CDR
[C, D, R] = cdr(F);
% Make D positive:
for k = 1:size(D)
    if ( D(k, k) < 0 )
        D(k, k) = -D(k, k);
        C(:, k) = -C(:,k);
    end
end

G1 = (C*D)*R.';
pass(3) = ( norm(G1 - F) < tol );
pass(4) = ( norm(F - G1) < tol );

G2 = C*(D*R.');
pass(5) = ( norm(G2 - F) < tol );
pass(6) = ( norm(F - G2) < tol );

G3 = (C*sqrt(D))*(sqrt(D)*R.');
pass(7) = ( norm(G3 - F) < tol );
pass(8) = ( norm(F - G3) < tol );

end


text:
function [A,b] = parallax(n)
%PARALLAX Stellar parallax problem with 28 fixed, real observations.
%
% [A,b] = parallax(n)
%
% Stellar parallax problem with 28 fixed, real observations.
%
% The underlying problem is a Fredholm integral equation of the
% first kind with kernel
%    K(s,t) = (1/sigma*sqrt(2*pi))*exp(-0.5*((s-t)/sigma)^2) ,
% and it is discretized by means of a Galerkin method with n
% orthonormal basis functions.  The right-hand side consists of
% a measured distribution function of stellar parallaxes, and its
% length is fixed, m = 26.  The exact solution, which represents
% the true distribution of stellar parallaxes, in not known.

% Reference: W. M. Smart, "Stellar Dynamics", Cambridge University Press,
% 1938; p. 30.

% Discretized by Galerkin method with orthonormal box functions;
% 2-D integration is done by means of the computational molecule:
%       1   4   1
%       4  16   4
%       1   4   1

% Per Christian Hansen, IMM, 09/16/92.

% Initialization.
a = 0; b = 0.1; m = 26; sigma = 0.014234;
hs = 0.130/m; hx = (b-a)/n; hsh = hs/2; hxh = hx/2;
ss = (-0.03 + [0:m-1]'*hs)*ones(1,n);
xx = ones(m,1)*(a + [0:n-1]*hx);

% Set up the matrix.
A =     16*exp(-0.5*((ss+hsh - xx-hxh)/sigma).^2);
A = A + 4*(exp(-0.5*((ss+hsh - xx    )/sigma).^2) + ...
           exp(-0.5*((ss+hsh - xx-hx )/sigma).^2) + ...
           exp(-0.5*((ss     - xx-hxh)/sigma).^2) + ...
           exp(-0.5*((ss+hs  - xx-hxh)/sigma).^2));
A = A +   (exp(-0.5*((ss     - xx    )/sigma).^2) + ...
           exp(-0.5*((ss+hs  - xx    )/sigma).^2) + ...
           exp(-0.5*((ss     - xx-hx )/sigma).^2) + ...
           exp(-0.5*((ss+hs  - xx-hx )/sigma).^2));
A = sqrt(hs*hx)/(36*sigma*sqrt(2*pi))*A;

% Set up the normalized right-hand side.
b = [3;7;7;17;27;39;46;51;56;50;43;45;43;32;33;29;...
     21;12;17;13;15;12;6;6;5;5]/(sqrt(hs)*640);

repo: REST


text:
f = [110*1.3 30*2.0 125*1.56 75*1.8 95*.95 100*2.25 50*1.35];
A = [120 210 150.75 115 186 140 85;
     110 30 125 75 95 100 50;
     1 1 1 1 1 1 1;
     1 -1 0 0 0 0 0;
     0 0 1 0 -2 0 0;
     0 0 0 -1 0 -1 1];
b = [55000;40000;400;0;0;0];
lp = lp_maker(f, A, b, [-1; -1; -1; -1; -1; -1], [10 10 10 10 20 20 20], [100 Inf 50 Inf Inf 250 Inf], [], 1, 0);
solvestat = mxlpsolve('solve', lp)
format bank
obj = mxlpsolve('get_objective', lp)
format short
x = mxlpsolve('get_variables', lp)
mxlpsolve('delete_lp', lp);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% Author: Jai Juneja
% Date: 22/03/2013
%
% Find the shortest Euclidean distance from a point to a line.
% Inputs:
%   pt  :   point in space
%   v1  :   first point along line
%   v2  :   second point along line
%           pt, v1 and v2 must all be at least 3x1 vectors
% Outputs:
%   d   :   shortest distance from point pt to line
function d = pointToLine(pt, v1, v2)
    a = v1 - v2;
    b = pt - v2;
    d = norm(cross(a,b)) / norm(a);
end

author: jaijuneja

repo: ekf-slam-matlab


text:
function M = spectrahedronfactory(n, k)
% Manifold of n-by-n symmetric positive semidefinite matrices of rank k
% with trace (sum of diagonal elements) equal to 1.
%
% function M = spectrahedronfactory(n, k)
%
% A point X on the manifold is parameterized as YY^T where Y is a matrix of
% size nxk. As such, X is symmetric, positive semidefinite. We restrict to
% full-rank Y's, such that X has rank exactly k. The point X is numerically
% represented by Y (this is more efficient than working with X, which may
% be big). Tangent vectors are represented as matrices of the same size as
% Y, call them Ydot, so that Xdot = Y Ydot' + Ydot Y and trace(Xdot) == 0.
% The metric is the canonical Euclidean metric on Y.
%
% The trace constraint on X (trace(X) == 1) translates to a unit Frobenius
% norm constraint on Y: trace(X) = norm(Y, 'fro')^2 == 1. The set of such
% Y's forms the unit sphere in R^(nxk): see spherefactory. But because for
% any orthogonal Q of size k, it holds that (YQ)(YQ)' = YY', we "group" all
% matrices of the form YQ in an equivalence class. The set of equivalence
% classes is a Riemannian quotient manifold, implemented here.
%
%
% Note that this geometry formally breaks down at rank-deficient Y's.
% As an alternative, you may use the sphere manifold (it has larger
% dimension (by 1), but does not break down at rank drop.)
%
% The geometry is taken from the 2010 paper:
% M. Journee, P.-A. Absil, F. Bach and R. Sepulchre,
% "Low-Rank Optimization on the Cone of Positive Semidefinite Matrices".
% Paper link: http://www.di.ens.fr/~fbach/journee2010_sdp.pdf
%
%
% Please cite the Manopt paper as well as the research paper:
%     @Article{journee2010low,
%       Title   = {Low-rank optimization on the cone of positive semidefinite matrices},
%       Author  = {Journ{\'e}e, M. and Bach, F. and Absil, P.-A. and Sepulchre, R.},
%       Journal = {SIAM Journal on Optimization},
%       Year    = {2010},
%       Number  = {5},
%       Pages   = {2327--2351},
%       Volume  = {20},
%       Doi     = {10.1137/080731359}
%     }
%
%
% See also: spherefactory elliptopefactory symfixedrankYYfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, July 11, 2013.
% Contributors: Nicolas Boumal
% Change log:
%
%   Apr. 2, 2015 (NB):
%       Replaced trace(A'*B) by A(:)'*B(:) (equivalent but faster).
%       Updated documentation.
%
%   Apr. 17, 2018 (NB):
%       Removed dependence on lyap.
%
%   Sep.  6, 2018 (NB):
%       Removed M.exp() as it was not implemented.

    M.name = @() sprintf('YY'' quotient manifold of %dx%d psd matrices of rank %d with trace 1', n, k);

    M.dim = @() n*k - 1 - k*(k-1)/2;

    % Euclidean metric on the total space
    M.inner = @(Y, eta, zeta) eta(:)'*zeta(:);

    M.norm = @(Y, eta) sqrt(M.inner(Y, eta, eta));

    M.dist = @(Y, Z) error('spectrahedronfactory.dist not implemented yet.');

    M.typicaldist = @() 10*k;

    M.proj = @projection;
    function etaproj = projection(Y, eta)
        % Projection onto the tangent space, i.e., on the tangent space of
        % ||Y|| = 1

        eta = eta - (eta(:)'*Y(:))*Y;

        % Projection onto the horizontal space
        YtY = Y'*Y;
        SS = YtY;
        AS = Y'*eta - eta'*Y;
        % Omega = lyap(SS, -AS);
        Omega = lyapunov_symmetric(SS, AS);
        etaproj = eta - Y*Omega;
    end

    M.tangent = M.proj;
    M.tangent2ambient = @(Y, eta) eta;

    M.retr = @retraction;
    function Ynew = retraction(Y, eta, t)
        if nargin < 3
            t = 1.0;
        end
        Ynew = Y + t*eta;
        Ynew = Ynew/norm(Ynew, 'fro');
    end

    M.egrad2rgrad = @(Y, eta) eta - (eta(:)'*Y(:))*Y;

    M.ehess2rhess = @ehess2rhess;
    function Hess = ehess2rhess(Y, egrad, ehess, eta)

        % Directional derivative of the Riemannian gradient
        Hess = ehess - (egrad(:)'*Y(:))*eta - ( (ehess(:)'*Y(:)) + (eta(:)'*egrad(:)) )*Y;
        Hess = Hess - (Hess(:)'*Y(:))*Y;

        % Project on the horizontal space
        Hess = M.proj(Y, Hess);

    end

    % Notice that the hash of two equivalent points will be different...
    M.hash = @(Y) ['z' hashmd5(Y(:))];

    M.rand = @random;

    function Y = random()
        Y = randn(n, k);
        Y = Y/norm(Y,'fro');
    end

    M.randvec = @randomvec;
    function eta = randomvec(Y)
        eta = randn(n, k);
        eta = projection(Y, eta);
        nrm = M.norm(Y, eta);
        eta = eta / nrm;
    end

    M.lincomb = @matrixlincomb;

    M.zerovec = @(Y) zeros(n, k);

    M.transp = @(Y1, Y2, d) projection(Y2, d);

    M.vec = @(Y, u_mat) u_mat(:);
    M.mat = @(Y, u_vec) reshape(u_vec, [n, k]);
    M.vecmatareisometries = @() true;

end

author: NicolasBoumal

repo: manopt


text:
%% build_system
% Build a system to solve by MaxwellFDFD.

%%% Syntax
%  [osc, grid3d, s_factor_cell, eps_cell, mu_cell, J_cell] = build_system(ge, OSC, DOM, OBJ, SRC, [progmark])
%  [..., obj_array, src_array, mat_array] = build_system(ge, OSC, DOM, OBJ, SRC, [pragmark])
%  [..., eps_node_array, mu_node_array] = build_system(ge, OSC, DOM, OBJ, SRC, [pragmark])

%%% Description
% |build_system(ge, OSC, DOM, OBJ, SRC, [progmark])| constructs a system from
% given objects and sources.  The constructed system is typically used inside
% <maxwell_run.html maxwell_run>.
%
% |ge| is an instance of |GT| and indicates the grid type of the _E_-field.
% Each following argument, |OSC|, |DOM|, |OBJ|, and |SRC|, represents a
% group of parameters. Each group supports several flexible expressions.
% For more details, see the relevant sections about the input parameter
% groups in <maxwell_run.html |maxwell_run|>.
%
% An additional input parameter |progmark| is an instance of <ProgMark.html
% ProgMark>, which outputs the progress of the system build procedure as the
% standard output.  If it is not given, then it is created internally.
%
% |[osc, grid3d, s_factor_cell, eps_cell, mu_cell, J_cell, M_cell] = build_system(...)|
% returns
%
% * |osc|, an instance of <Oscillation.html Oscillation>
% * |grid3d|, an instance of <Grid3d.html Grid3d>,
% * |s_factor_cell|, a cell array of PML s-factors: |{sx_array, sy_array,
% sz_array}|
% * |eps_cell|, a cell array of electric permittivity evaluated at the E-field
% positions: |{eps_xx_array, eps_yy_array, eps_zz_array}|
% * |mu_cell|,  a cell array of magnetic permeability evaluated at the H-field
% positions: |{mu_xx_array, mu_yy_array, mu_zz_array}|
% * |J_cell|, a cell array of electric current sources: |{Jx_array, Jy_array,
% Jz_array}|
% * |M_cell|, a cell array of electric current sources: |{Mx_array, My_array,
% Mz_array}|
%
% |[..., obj_array, src_array, mat_array] = build_system(...)| returns
% additionally arrays of instances of <EMObject.html |EMObject|>, <Source.html
% |Source|>, and <Material.html |Material|>.  The |EMObject| and |Source|
% elements represent the objects and sources placed in the simulation
% domain, so they can be used to visualize the simulation domain.
%
% |[..., eps_node_array, mu_node_array] = build_system(...)| returns
% additionally arrays of electric permittivity and magnetic permeability
% evaluated at the nodes of the finite-difference grid.

%%% Example
%   gray = [0.5 0.5 0.5];  % [r g b]
%   [osc, grid3d, s_factor_cell, eps_cell, mu_cell, J_cell, M_cell,	...
%       obj_array, src_array, mat_array, eps_node_array, mu_node_array] = build_system(...
%       'OSC', 1e-9, 1550, ...
%       'DOM', {['Palik/SiO2'], 'none'}, [-700, 700; -600, 600; -200, 1700], 20, BC.p, 200, ...
%       'OBJ', ...
%           {['Palik/SiO2'], 'none'}, Box([-50, 50; -50, 50; -200, 1700], [2, 2, 20]), ...  % OBJ1
%           {['CRC/Ag'], gray}, [Box([-700, -25; -25, 25; -200, 1700], 20), Box([25, 700; -25, 25; -200, 1700], 20)], ...  % OBJ2
%       'SRC', PointSrc(Axis.x, [0, 0, 200]) ...
%       );

function [osc, grid3d, s_factor_cell, eps_cell, mu_cell, J_cell, M_cell, ...
	obj_array, src_array, mat_array, eps_node, mu_node, isiso] = build_system(varargin)

	iarg = nargin; arg = varargin{iarg};
	if istypesizeof(arg, 'ProgMark')
		pm = arg;
		narglim = nargin - 1;
	else
	 	pm = ProgMark();
		varargin = [varargin, {pm}];
		narglim = nargin;
	end

	iarg = 1; arg = varargin{iarg};
	chkarg(istypesizeof(arg, 'GT'), 'argument #%d should be "ge" (GT).', iarg);
	ge = arg;

	iarg = iarg + 1; arg = varargin{iarg};
	chkarg(istypesizeof(arg, 'PML'), 'argument #%d should be "pml" (PML).', iarg);
	pml = arg;

	function material = create_material(varargin)
		narg = nargin;
		if istypesizeof(varargin{end}, 'logical')
			narg = narg - 1;
		end

		chkarg(narg >= 2, '# of arguments should be at least 2.')
		matname = varargin{1};
		if isempty(strfind(matname, '/'))  % data table is not specified
			material = Material(varargin{:});
		else  % data table is specified
			material = Material.fromtable(osc, varargin{:});
		end
	end

	osc = Oscillation.empty();
	obj_dom = EMObject.empty();
	shape_array = Shape.empty();
	sshape_array = Shape.empty();
	mat_array = Material.empty();
	obj_array = EMObject.empty();
	sobj_array = EMObject.empty();
	srcj_array = [];
	srcm_array = [];
	isepsgiven = false;
	isTFSF = false;
	while iarg < narglim
		iarg = iarg + 1; arg = varargin{iarg};

		if ischar(arg) && strcmpi(arg,'OSC')
			% Set up OSC.
			iarg = iarg + 1; arg = varargin{iarg};
			chkarg((istypesizeof(arg, 'real') && arg > 0) || istypesizeof(arg, 'Oscillation'), ...
				'"argument #%d should be either "L0" (positive) or "osc" (instance of Oscillation).', iarg);
			if istypesizeof(arg, 'real')
				L0 = arg;
				iarg = iarg + 1; wvlen = varargin{iarg};
				chkarg(istypesizeof(wvlen, 'complex'), 'argument #%d should be "wvlen" (complex).', iarg);
				unit = PhysUnit(L0);
				osc = Oscillation(wvlen, unit);
			else  % arg is instance of Oscillation
				osc = arg;
			end
		elseif ischar(arg) && strcmpi(arg,'DOM')
			% Set up DOM.
			iarg = iarg + 1; arg = varargin{iarg};
			chkarg(iscell(arg) || istypesizeof(arg, 'Material') || istypesizeof(arg, 'EMObject'), ...
				'argument #%d should be cell, instance of Material, or instance of EMObject.', iarg);
			if istypesizeof(arg, 'EMObject')
				obj_dom = arg;
				domain = obj_dom.shape;
				chkarg(istypesizeof(domain, 'Domain'), 'argument #%d should be instance of EMObject with Domain as its shape.', iarg);
			else
				if iscell(arg)
					mat_dom = create_material(arg{:});
				else
					assert(istypesizeof(arg, 'Material'));
					mat_dom = arg;
				end

				iarg = iarg + 1; arg = varargin{iarg};
				chkarg(istypesizeof(arg, 'real', [Axis.count, Sign.count]) || istypesizeof(arg, 'Domain'), ...
					'argument #%d should be either "box_dom" ([xmin xmax; ymin ymax; zmin zmax]) or "domain" (instance of Domain).', iarg);
				if istypesizeof(arg, 'real', [Axis.count, Sign.count])
					box_domain = arg;
					iarg = iarg + 1; dl_domain = varargin{iarg};
					chkarg(istypeof(dl_domain, 'real') && isexpandable2row(dl_domain, Axis.count), ...
						'"argument #%d should be dl_domain (positive number or length-%d row vector of positive numbers).', iarg, Axis.count);
					domain = Domain(box_domain, expand2row(dl_domain, Axis.count));
				else  % arg is instance of Domain
					domain = arg;
				end
				obj_dom = EMObject(domain, mat_dom);
			end
			mat_array = [mat_array(1:end), obj_dom.material];

			% Set up boundary conditions and PML thicknesses.
			iarg = iarg + 1; bc = varargin{iarg};
			chkarg(istypeof(bc, 'BC') && isexpandable2mat(bc, Axis.count, Sign.count), ...
				'argument #%d should be "bc" (scalar, length-%d row vector, or %d-by-%d matrix with BC as elements).', iarg, Axis.count, Axis.count, Sign.count);
			iarg = iarg + 1; Lpml = varargin{iarg};
			chkarg(istypeof(Lpml, 'real') && isexpandable2mat(Lpml, Axis.count, Sign.count) && all(all(Lpml>=0)), ...
				'argument #%d should be "Lpml" (scalar, length-%d row vector, or %d-by-%d matrix with nonnegative numbers as elements).', iarg, Axis.count, Axis.count, Sign.count);

			% Set up the degree of the polynomial grading of the PML scale
			% factors.
			iarg = iarg + 1; arg = varargin{iarg};
			deg_pml = 4;  % polynomial degree
			if istypeof(arg, 'real')
				deg_pml = arg;
			else
				iarg = iarg - 1; % because deg_pml is optional argument
			end

			% Set up the target reflection coefficient of the PML.
			iarg = iarg + 1; arg = varargin{iarg};
			R_pml = exp(-16);  % target reflection coefficient
			if istypeof(arg, 'real')
				R_pml = arg;
			else
				iarg = iarg - 1; % because R_pml is optional argument
			end

			% Set up a flag to generate a grid dynamically.
			iarg = iarg + 1; arg = varargin{iarg};
			withuniformgrid = false;  % generate a grid dynamically by default.
			if istypesizeof(arg, 'logical')
				withuniformgrid = arg;
			else
				iarg = iarg - 1; % because withuniformgrid is optional argument
			end
		elseif ischar(arg) && (strcmpi(arg,'OBJ') || strcmpi(arg,'SOBJ'))
			% Set up OBJ.
			is_scatterer = strcmpi(arg,'SOBJ');
			iarg = iarg + 1; arg = varargin{iarg};
			if istypesizeof(arg, 'complex', [0 0 0])  % 3D complex array with arbitrary size
				isepsgiven = true;
				eps_node_cell = {arg, arg, arg};
				mu_node_temp = ones(size(arg));
				mu_node_cell = {mu_node_temp, mu_node_temp, mu_node_temp};
			else
				% Set up objects.
				obj_array_temp = EMObject.empty();
				shape_array_temp = Shape.empty();
				while iscell(arg) || istypesizeof(arg, 'Material') || istypesizeof(arg, 'EMObject', [1 0])
					if istypesizeof(arg, 'EMObject', [1 0])
						objs = arg;
						obj_array_temp = [obj_array_temp(1:end), objs];
						for obj = objs
							shape_array_temp = [shape_array_temp(1:end), obj.shape];
							mat_array = [mat_array(1:end), obj.material];
						end
						iarg = iarg + 1; arg = varargin{iarg};
					else
						if iscell(arg)
							mat = create_material(arg{:});
						else
							assert(istypesizeof(arg, 'Material'));
							mat = arg;
						end
						mat_array = [mat_array(1:end), mat];

						iarg = iarg + 1; arg = varargin{iarg};
						while istypesizeof(arg, 'Shape', [1 0])
							shapes = arg;
							shape_array_temp = [shape_array_temp(1:end), shapes];
							objs = EMObject(shapes, mat);
							obj_array_temp = [obj_array_temp(1:end), objs];
							iarg = iarg + 1; arg = varargin{iarg};
						end
					end
				end
				iarg = iarg - 1;

				if is_scatterer
					sshape_array = [sshape_array(1:end), shape_array_temp];
					sobj_array = [sobj_array(1:end), obj_array_temp];
				else
					shape_array = [shape_array(1:end), shape_array_temp];
					obj_array = [obj_array(1:end), obj_array_temp];
				end
			end
		elseif ischar(arg) && strcmpi(arg,'SRCJ')
			% Set up sources.
			iarg = iarg + 1; arg = varargin{iarg};
			if ~istypesizeof(arg, 'Source', [1 0])
				warning('Maxwell:buildSys', 'no source is given.');
			end

			while istypesizeof(arg, 'Source', [1 0])
				if istypesizeof(arg, 'TFSFPlaneSrc')
					isTFSF = true;
				end
				srcj_array_curr = arg;
				for src = srcj_array_curr
					src.set_gridtype(ge);
				end
				srcj_array = [srcj_array(1:end), srcj_array_curr];
				iarg = iarg + 1; arg = varargin{iarg};
			end
			iarg = iarg - 1;
		elseif ischar(arg) && strcmpi(arg,'SRCM')
			% Set up sources.
			iarg = iarg + 1; arg = varargin{iarg};
			if ~istypesizeof(arg, 'Source', [1 0])
				warning('Maxwell:buildSys', 'no source is given.');
			end

			while istypesizeof(arg, 'Source', [1 0])
				if istypesizeof(arg, 'TFSFPlaneSrc')
					isTFSF = true;
				end
				srcm_array_curr = arg;
				for src = srcm_array_curr
					src.set_gridtype(alter(ge));
				end
				srcm_array = [srcm_array(1:end), srcm_array_curr];
				iarg = iarg + 1; arg = varargin{iarg};
			end
			iarg = iarg - 1;
		elseif iarg == narglim
			chkarg(false, ['some arguments are not used.\n', ...
				'Suggestion: check if each parameter group is specified with beginning specifier.']);
		end
	end
	chkarg(~isempty(osc), 'OSC parameter groups should be set.');
	chkarg(~isempty(obj_dom), 'DOM parameter groups should be set.');
	obj_array = [obj_dom, obj_array];
	src_array = [srcj_array, srcm_array];
	if isTFSF && isempty(sobj_array)
		warning('Maxwell:objAssign', 'TF/SF source is used, but scatteres are not defined in SOBJ group.');
	end

	chkarg(iarg <= narglim, 'more arguments than expected.');
	pm.mark('initial setup');

	fprintf('\tLength Unit: %s m\n', osc.unit.value(PhysQ.L));
	fprintf('\twvlen = %s, freq = %s eV\n', num2str(osc.in_L0()), num2str(osc.in_eV()));

	mat_array = unique(mat_array);
	isiso = true;
	fprintf('materials used:\n');
	for mat = mat_array
		epstext = mat.eps;
		if length(unique(epstext)) == 1
			epstext = epstext(Axis.x);
		end

		mutext = mat.mu;
		if length(unique(mutext)) == 1
			mutext = mutext(Axis.x);
		end

		fprintf('\t%s: eps = %s, mu = %s\n', mat.name, mat2str(epstext), mat2str(mutext));

		isiso = isiso && mat.isiso;
	end

	% Generate a grid.
	[lprim, Npml] = generate_lprim3d(domain, Lpml, [shape_array, sshape_array], src_array, withuniformgrid);
	grid3d = Grid3d(osc.unit, lprim, Npml, bc);
	if withuniformgrid
		pm.mark('uniform grid generation');
	else
		pm.mark('nonuniform grid generation');
	end
	fprintf('\t[Nx Ny Nz] = %s\n', mat2str(grid3d.N));

	% Generate a warning when a seemingly 2D simulation is defined on a 3D grid.
	[like2d, normal_axis] = is2dlike(grid3d.N);
	if like2d && grid3d.N(normal_axis) >= 2  % possible user mistake
		warning(['If this a 2D structure, N%s should be 1; ', ...
			'check d%s''s of objects and locations of sources'], char(normal_axis), char(normal_axis));
	end

	% Construct material parameters.
	if ~isepsgiven
		[eps_node_cell, mu_node_cell] = assign_material_node(grid3d, obj_array);  % Nx x Ny x Nz
	end
	eps_cell = mean_material_node(grid3d, ge, eps_node_cell);
	mu_cell = mean_material_node(grid3d, alter(ge), mu_node_cell);

	% Construct PML s-factors.
	s_factor_cell = generate_s_factor(osc.in_omega0(), grid3d, deg_pml, R_pml);
	pm.mark('eps and mu assignment');

	if ~isTFSF
		% Solve for modes.
		for src = src_array
			if istypesizeof(src, 'ModalSrc')
				modalsrc = src;
				if ~modalsrc.ispreped
					prep_modalsrc(ge, pml, osc, grid3d, eps_cell, mu_cell, s_factor_cell, modalsrc);
				end

				neff = modalsrc.neff;
				beta = 2*pi*neff / osc.in_L0();
				pm.mark('mode calculation');
				fprintf('\tbeta = %s, n_eff = %s\n', num2str(beta), num2str(neff));
			end
		end
	else  % isTFSF == true
		% Set up J for TF/SF.
		for src = src_array
			if istypesizeof(src, 'TFSFPlaneSrc')
				tfsfsrc = src;
				cb_center = num2cell(tfsfsrc.shape.cb_center);
				for bgobj = fliplr(obj_array)
					if bgobj.shape.contains(cb_center{:})
						break;  % assume that last object containing TF box center fills TF box
					end
				end
				tfsfsrc.set_bg_material(bgobj.material);
				F0 = tfsfsrc.create_incidentF(osc, grid3d);
				JM = cell(1, Axis.count);
				for w = Axis.elems
					JM{w} = zeros(grid3d.N);
				end

				if tfsfsrc.gt == ge
					eqtype_tfsf = EquationType(FT.e, ge);  % for SRCJ, create E-field eq
				else
					eqtype_tfsf = EquationType(FT.h, ge);  % for SRCM, create H-field eq
				end
% 				A = create_eq(eqtype_tfsf, pml, osc.in_omega0(), eps_cell, mu_cell, s_factor_cell, JM, JM, grid3d);
				eq = MatrixEquation(eqtype_tfsf, pml, osc.in_omega0(), eps_cell, mu_cell, s_factor_cell, JM, JM, grid3d);
				Op = eq.matrixfree_op();

				x0 = [F0{Axis.x}(:); F0{Axis.y}(:); F0{Axis.z}(:)];
				r = reordering_indices(Axis.count, grid3d.N);
				x0 = x0(r);

% 				JM = (A*x0) ./ (-1i*osc.in_omega0());
				JM = Op(x0, 'notransp') ./ (-1i*osc.in_omega0());
				JM = reshape(JM, [Axis.count grid3d.N]);
				JM = permute(JM, [Axis.elems+1, 1]);
				JM = {JM(:,:,:,Axis.x), JM(:,:,:,Axis.y), JM(:,:,:,Axis.z)};

				tfsfsrc.setJM(JM, grid3d);
			end
		end

		% Add sobj_array to the already-generated eps and mu.
		[eps_node_cell, mu_node_cell] = assign_material_node(grid3d, sobj_array, eps_node_cell, mu_node_cell);  % Nx x Ny x Nz
		eps_cell = mean_material_node(grid3d, ge, eps_node_cell);  % Nx x Ny x Nz
		mu_cell = mean_material_node(grid3d, alter(ge), mu_node_cell);  % Nx x Ny x Nz

		pm.mark('TF/SF source assignment');
	end
	obj_array = [obj_array, sobj_array];

	eps_node = cell(1, Axis.count);
	mu_node = cell(1, Axis.count);
	for w = Axis.elems
		eps_node_cell{w} = expand_node_array(grid3d, eps_node_cell{w});  % (Nx+2) x (Ny+2) x (Nz+2)
		mu_node_cell{w} = expand_node_array(grid3d, mu_node_cell{w});  % (Nx+2) x (Ny+2) x (Nz+2)

		eps_node{w} = Scalar3d(eps_node_cell{w}, grid3d, [GT.dual GT.dual GT.dual], osc, PhysQ.eps, '\epsilon');
		mu_node{w} = Scalar3d(mu_node_cell{w}, grid3d, [GT.dual GT.dual GT.dual], osc, PhysQ.mu, '\mu');
	end

	% Construct sources.
	J_cell = assign_source(grid3d, srcj_array);
	M_cell = assign_source(grid3d, srcm_array);
	pm.mark('J assignment');
end


text:
function neighborsCell = identifyNeighborsAuto(inputImages, inputSignals, varargin)
	% This code automatically sorts through to find all obj neighbors within a certain distance of the target (boundary to boundary). The output is a cell array with vectors of the neighbor indices to each obj.
	% Biafra Ahanonu
	% started: 2013.11.01
	% based on code by laurie burns, started: sept 2010.
	% inputs
		% inputImages - [x y nCells] matrices containing each set of filters
		% inputSignals - [nFilters frames] matrices containing each set of filter traces
	% options
		% _
	% outputs
		% _

	[neighborsCell] = ciapkg.neighbor.identifyNeighborsAuto(inputImages, inputSignals,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah