text:
% Pattern Recognition Tools (PRTools Guide)
% Version 5.1.1 14-May-2014
%
%Datasets and Mappings (just most important routines)
%---------------------
%prdataset      - Define dataset from datamatrix and labels
%datasets       - List information on datasets (just help, no command)
%prdatafile     - Define dataset from directory of object files
%datafiles      - List information on datafiles (just help, no command)
%cat2data       - Create categorical dataset
%classnames     - Retrieve names of classes
%classsizes     - Retrieve sizes of classes
%feat2lab       - Label dataset by one of its features and remove this feature
%gencirc        - Generation of a one-class circular dataset
%genclass       - Generate class frequency distribution
%genlab         - Generate dataset labels
%getlab         - Retrieve object labels from datasets and mappings
%getnlab        - Retrieve nummeric object labels from dataset
%setfeatlab     - Set feature labels in dataset
%getfeatlab     - Get feature labels in dataset
%getfeat        - Retrieve feature labels from datasets and mappings
%setdat         - Change data in dataset for classifier output
%setdata        - Change data in dataset or mapping
%getdata        - Retrieve data from dataset or mapping
%setlabels      - Change labels of dataset or mapping
%getlabels      - Retrieve labels from a dataset
%setprior       - Reset class prior probabilities of dataset
%getprior       - Retrieve class prior probabilities from dataset
%addlabels      - Add additional labelling
%changelablist  - Change current active labeling
%misval         - Fix missing values in a dataset
%multi_labeling - List information on multi-labeling (help only)
%prmapping      - Define and retrieve mapping and classifier from data
%mappings       - List information on mappings (just help, no command)
%renumlab       - Convert labels to numbers
%matchlab       - Match different labelings
%prarff         - Convert ARFF file (WEKA) to PRTools dataset
%remclass       - Remove a class from a dataset
%seldat         - Retrieve a part of a dataset
%selclass       - Retrieve a class from a dataset
%
%Data Generation (more in prdatasets)
%---------------
%circles3d   - Create a dataset containing 2 circles in 3 dimensions
%lines5d     - Create a dataset containing 3 lines in 5 dimensions
%gendat      - Random sampling of datasets for training and testing
%gensubsets  - Generation of a consistent series of subsets of a dataset
%gendatgauss - Generation of multivariate Gaussian distributed data
%gendatb     - Generation of banana shaped classes
%gendatc     - Generation of circular classes
%gendatd     - Generation of two difficult classes
%gendath     - Generation of Highleyman classes
%gendati     - Generation of random windows from images
%gendatk     - Nearest neighbour data generation
%gendatl     - Generation of Lithuanian classes
%gendatm     - Generation of 8 2d classes
%gendatp     - Parzen density data generation
%gendatr     - Generate regression dataset from data and target values
%gendats     - Generation of two Gaussian distributed classes
%gendatw     - Sample dataset by given weigths
%gendatv     - Generation of a very large dataset
%gentrunk    - Generation of Trunk's example
%prdata      - Read data from file
%seldat      - Select classes / features / objects from dataset
%spirals     - Generation of a two-class spiral dataset
%getwindows  - Get pixel feature vectors around given pixels in image dataset
%prdataset   - Read existing dataset from file
%prdatasets  - Overview and download of standard datasets
%
%Datafiles
%---------
%prdatafile     - Define datafile from set of files in directory
%createdatafile - Save datafile, store intermediate result as raw datafile
%savedatafile   - Save datafile, store intermediate result as mature datafile
%filtm          - Mapping for arbitrary processing of a datafile
%prdatafiles    - Overview and download of standard datafiles
%
%Linear and Quadratic Classifiers (*operate on datasets and datafiles)
%--------------------------------
%fisherc     - Minimum least square linear classifier
%ldc         - Normal densities based linear (muli-class) classifier
%loglc       - Logistic linear classifier
%nmc         - Nearest mean linear classifier
%nmsc        - Scaled nearest mean linear classifier
%quadrc      - Quadratic classifier
%qdc         - Normal densities based quadratic (multi-class) classifier
%udc         - Uncorrelated normal densities based quadratic classifier
%klldc       - Linear classifier based on KL expansion of common cov matrix
%pcldc       - Linear classifier based on PCA expansion on the joint data
%polyc       - Add polynomial features and run arbitrary classifier
%subsc       - Subspace classifier
%statslinc   - Linear classifier from the Stats toolbox
%
%classc      - Converts a mapping into a classifier
%labeld      - Find labels of objects by classification
%logdens     - Convert density estimates to log-densities for more accuracy
%rejectc     - Creates reject version of exisiting classifier
%testc       - General error estimation routine for trained classifiers
%
%Other Classifiers
%-----------------
%knnc        - k-nearest neighbour classifier (find k, build classifier)
%testk       - Error estimation for k-nearest neighbour rule
%edicon      - Edit and condense training sets
%statsknnc   - k-nearest neighbour classifier from the Stats toolbox
%
%weakc       - Weak classifier
%stumpc      - Decision stump classifier
%adaboostc   - ADABoost classifier
%
%parzenc     - Parzen classifier
%parzendc    - Parzen density based classifier
%testp       - Error estimation for Parzen classifier
%
%treec       - Construct binary decision tree classifier
%dtc         - Decision tree classifier, rewritten, also for nominal features
%statsdtc    - Decision tree classifier from the Stats toolbox
%randomforestc - Breiman's random forest classifier
%naivebc     - Naive Bayes classifier
%statsnbc    - Naive Bayes classifier from the Stats toolbox
%bpxnc       - Feed forward neural network classifier by backpropagation
%lmnc        - Feed forward neural network by Levenberg-Marquardt rule
%neurc       - Automatic neural network classifier
%perlc       - Linear perceptron
%rbnc        - Radial basis neural network classifier
%rnnc        - Random neural network classifier
%ffnc        - Feed-forward neural net classifier back-end routine
%bagc        - Feature set classifier, e.g. for multiple-instance learning
%
%fdsc        - Feature based dissimilarity space classifier
%mdsc        - Manhatten distance feature based dissimilarity space classifier
%vpc         - Voted perceptron classifier
%drbmc       - Discriminative restricted Boltzmann machine classifier
%
%libsvc      - Support vector classifier by LIBSVM
%nulibsvc    - Support vector classifier by LIBSVM
%svc         - Support vector classifier
%svo         - Support vector optimizer
%nusvc       - Support vector classifier
%nusvo       - Support vector optimizer
%rbsvc       - Radial basis SV classifier
%kernelc     - General kernel/dissimilarity based classification
%
%Normal Density Based Classification
%-----------------------------------
%distmaha    - Mahalanobis distance
%meancov     - Estimation of means and covariance matrices from multiclass data
%nbayesc     - Bayes classifier for given normal densities
%ldc         - Normal densities based linear (muli-class) classifier
%qdc         - Normal densities based quadratic (multi-class) classifier
%udc         - Uncorrelated normal densities based quadratic classifier
%mogc        - Mixture of gaussians classification
%testn       - Error estimate of discriminant on normal distributions
%
%Feature Selection
%-----------------
%feateval    - Evaluation of a feature set
%featrank    - Ranking of individual feature permormances
%featsel     - Feature Selection
%featselb    - Backward feature selection
%featself    - Forward feature selection
%featsellr   - Plus-l-takeaway-r feature selection
%featseli    - Feature selection on individual performance
%featselm    - Feature selection map, general routine for feature selection
%featselo    - Branch and bound feature selection
%featselp    - Floating forward feature selection
%featselv    - Selection of varying features
%
%Classifiers and tests (general)
%-------------------------------
%bayesc      - Bayes classifier by combining density estimates
%classim     - Classify image using a given classifier
%classc      - Convert mapping to classifier
%labeld      - Find labels of objects by classification
%cleval      - Classifier evaluation (learning curve)
%clevalb     - Classifier evaluation (learning curve), bootstrap version
%clevalf     - Classifier evaluation (feature size curve)
%clevals     - Classifier evaluation (feature /learning curve), bootstrap
%confmat     - Computation of confusion matrix
%costm       - Cost mapping, classification using costs
%prcrossval  - Crossvalidation
%cnormc      - Normalisation of classifiers
%disperror   - Display error matrix with information on classifiers and datasets
%labelim     - Construct image of labeled pixels
%logdens     - Convert density estimates to log-densities for more accuracy
%loso        - Leave_one_set_out crossvalidation
%mclassc     - Computation of multi-class classifier from 2-class discriminants
%regoptc     - Optimisation of regularisation and complexity parameters
%reject      - Compute error-reject trade-off curve
%prroc       - Receiver-operator curve (ROC)
%shiftop     - Shift operating point of classifier
%testc       - General error estimation routine for trained classifiers
%testd       - Error of dataset applied to given classifier
%testauc     - Estimate error as area under the ROC
%
%Mappings
%--------
%affine      - Construct affine (linear) mapping from parameters
%bhatm       - Two-class Bhattacharryya mapping
%cmapm       - Compute some special maps
%datasetm    - Mapping conversion dataset
%disnorm     - Normalization of a dissimilarity matrix
%featselm    - Feature selection map, general routine for feature selection
%fisherm     - Fisher mapping
%chernoffm   - Chernoff mapping
%invsigm     - Inverse sigmoid map
%filtm       - Arbitrary operation on datafiles/datasets, object by object
%mapm        - Arbitrary mapping operation on doubles and datasets
%gaussm      - Mixture of Gaussians density estimation
%kernelm     - Kernel mapping
%klm         - Decorrelation and Karhunen Loeve mapping (PCA)
%klms        - Scaled version of klm, useful for prewhitening
%knnm        - k-Nearest neighbor density estimation
%mclassm     - Computation of mapping from multi-class dataset
%prmap       - General routine for computing and executing mappings
%mappingtools - Macro defining some mappings
%nlfisherm   - Nonlinear Fisher mapping
%normm       - Object normalization map
%parzenm     - Parzen density estimation
%parzenml    - Optimization of smoothing parameter in Parzen density estimation.
%pcam        - Principal Component Analysis
%pcaklm      - Backend routine for PC and KL mappings
%proxm       - Proximity mapping and kernel construction
%reducm      - Reduce to minimal space mapping
%remoutl     - Remove outliers
%rejectm     - Creates rejecting mapping
%scalem      - Compute scaling data
%sigm        - Simoid mapping
%spatm       - Augment image dataset with spatial label information
%tsnem       - tSNE mapping
%sammonm     - Multi-dimensional scaling by Sammon mapping
%userkernel  - User supplied kernel definition
%
%gtm         - Fit a Generative Topographic Mapping (GTM) by EM
%plotgtm     - Plot a Generative Topographic Mapping in 2D
%som         - Simple routine computing a Self-Organizing Map (SOM)
%prplotsom   - Plot a Self-Organizing Map in 2D
%
%Classifier combiners
%--------------------
%averagec    - Combining linear classifiers by averaging coefficients
%baggingc    - Bootstrapping and aggregation of classifiers
%dcsc        - Dynamic Classifier Selecting Combiner
%modselc     - Model Selection Combiner (Static selection)
%rsscc       - Random subspace combining classifier
%votec       - Voting classifier combiner
%wvotec      - Weighted voting classifier combiner
%maxc        - Maximum classifier combiner
%minc        - Minimum classifier combiner
%meanc       - Mean classifier combiner
%medianc     - Median classifier combiner
%mlrc        - Muli-response linear regression combiner
%naivebcc    - Naive Bayes classifier combiner
%perc        - Percentile combiner
%prodc       - Product classifier combiner
%traincc     - Train combining classifier
%fixedcc     - Fixed combiner construction, back end
%parsc       - Parse classifier or map
%rejectc     - Creates reject version of exisiting classifier
%parallel    - Parallel combining of classifiers
%bagcc       - Feature set combining classifier
%stacked     - Stacked combining of classifiers
%sequential  - Sequential combining of classifiers
%
%
%Regression
%----------
%linearr     - Linear regression
%ridger      - Ridge regression
%lassor      - LASSO
%svmr        - Support vector regression
%ksmoothr    - Kernel smoother
%knnr        - k-nearest neighbor regression
%pinvr       - Pseudo-inverse regression
%plsr        - Partial least squares regression
%plsm        - Partial least squares mapping
%gpr         - Gaussian Process regression
%
%testr       - Mean squared regression error
%rsquared    - R^2-statistic
%
%Handling images in datasets and datafiles
%-----------------------------------------
%data2im     - Convert dataset to image
%getobjsize  - Retrieve image size of feature images in datasets
%getfeatsize - Retrieve image size of object images in datasets
%obj2feat    - Transform object images to feature images in dataset
%feat2obj    - Transform feature images to object images in dataset
%im2feat     - Convert image to feature in dataset
%im2obj      - Convert image to object in dataset
%imsize      - Retrieve size of specific image in datafile
%im_patch    - Find / generate patches in object images
%band2obj    - Convert image bands to objects in dataset
%bandsel     - Select image bands in dataset or datafile
%selectim    - Select image in multi-band object image dataset/datafile
%show        - Display objects in datasets, datafiles and mappings
%im_dbr      - Image Database Retrieval GUI
%
%Operations on images in datasets and datafiles
%----------------------------------------------
%classim     - Classify image using a given classifier
%doublem     - Convert datafile images into double
%filtim      - Image operation on objects in datafiles/datasets
%spatm       - Augment image dataset with spatial label information
%im_box            - Bounding box
%im_center         - Center image
%im_fft            - FFT transform (and more)
%im_gauss          - Gaussian filtering by Matlab
%im_gray           - Multi-band to gray-value conversion
%im_hist_equalize  - Histogram equalization
%im_invert         - Invert image
%im_label          - Labeling binary images
%im_norm           - Normalize images w.r.t. mean and variance
%im_resize         - Resize images
%im_rotate         - Rotate images
%im_scale          - Scale images
%im_select_blob    - Select largest blob
%im_stretch        - Contrast stretching of images
%im_threshold      - Threshold images
%im_unif           - Uniform filtering
%
%Feature extraction from images in datasets and datafiles
%--------------------------------------------------------
%histm         - Convert images to histograms. Trains the bin positions
%im_hist       - Convert images to histograms for fixed bin positions
%im_harris     - Find Harris points in images
%im_moments    - Computes moments as features from object images
%im_mean       - Computes center of gravity
%im_measure    - Computes some measurements
%im_profile    - Computes image profiles
%im_skel_meas  - Skeleton measurements
%im_stat       - Compute some simple statistics
%
%Clustering and distances
%------------------------
%distm       - Distance matrix between two data sets
%emclust     - Expectation - maximization clustering
%proxm       - Proximity mapping and kernel construction
%hclust      - Hierarchical clustering
%kcentres    - k-centres clustering
%prkmeans    - k-means clustering
%modeseek    - Clustering by modeseeking
%
%mds         - Non-linear mapping by multi-dimensional scaling (Sammon)
%mds_cs      - Linear mapping by classical scaling
%mds_init    - Initialisation of multi-dimensional scaling
%mds_stress  - Dissimilarity of distance matrices
%
%Plotting
%--------
%gridsize    - Set gridsize used in the PRTools plot commands
%plotc       - Plot discriminant function for two features
%plote       - Plot error curves
%plotf       - Plot feature distribution
%plotm       - Plot mapping
%ploto       - Plot object functions
%plotr       - Plot regression functions
%plotdg      - Plot dendrgram (see hclust)
%scatterd    - Scatterplot
%scatterdui  - Scatterplot scatterplot with feature selection
%scattern    - Simple, unannotated scatterplot, no axes.
%scatterr    - Scatter regression dataset
%
%Various tests and support routines
%----------------------------------
%cdats              - Support routine for checking datasets
%concatm            - Concatenate cell array of mappings or datasets ({} --> [])
%iscomdset          - Test on compatible datasets
%isdataim           - Test on image dataset
%isdataset          - Test on dataset
%isfeatim           - Test on feature image dataset
%ismapping          - Test on mapping
%isobjim            - Test on object image dataset
%issequential       - Test on sequential mapping
%isstacked          - Test on stacked mapping
%isparallel         - Test on parallel mapping
%issym              - Test on symmetric matrix
%isvaldset          - Test on valid dataset
%isvaldfile         - Test on valid datafile
%matchlablist       - Match entries of label lists
%mapex              - Train and execute mapping on the same dataset
%labcmp             - Compare two label lists and find the differences
%nlabcmp            - Compare two label lists and count the differences
%testdatasize       - Check datasize and convert datafile to dataset
%define_mapping     - Define empty mapping
%mapping_task       - Check mapping task
%trained_mapping    - Defined trained mapping
%trained_classifier - Define trained classifier
%setdefaults        - Substitute defaults
%shiftargin         - Conditional shift of input arguments
%prload             - Load prtools4 mat-files and convert to prtools5
%prtools4to5        - Convert prtools4 directory to prtools5
%
%Examples
%--------
%prex_cleval     - learning curves
%prex_combining  - classifier combining
%prex_confmat    - confusion matrix, scatterplot and gridsize
%prex_datafile   - datafile usage
%prex_datasets   - standard datasets
%prex_density    - Various density plots
%prex_eigenfaces - Use of images and eigenfaces
%prex_matchlab   - K-means clustering and matching labels
%prex_mcplot     - Multi-class classifier plot
%prex_plotc      - Dataset scatter and classifier plot
%prex_mds        - Multi-dimensional scaling and visualisation
%prex_som        - Training a SelfOrganizing Maps
%prex_spatm      - Spatial smoothing of image classification
%prex_cost       - Cost matrices and rejection
%prex_logdens    - Density based classifier improvement
%prex_soft       - Soft label example
%prex_regr       - Regression example
%
%prdownload  - low level routine for retrieving datasets
%prglobal    - set / list all globals and settings
%prversion   - returns version information on PRTools
%prwaitbar   - report PRTools progress by single waitbar
%prwarning   - control PRTools warning level
%prmemory    - controol PRTools large dataset handling
%prtver      - prtools version back end
%typp        - list prtools routine nicely
%
%--- PRTools Guide ---

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function Cy = computeCy(ny,n_state)
    Cy_tmp1 = cell(ny,ny);
    Cy_tmp2 = cell(ny,ny );
    %Avant 5 a la place de 7
    for i=1:ny
       for j=1:ny
          if i==j
              Cy_tmp1{i,j} = eye(n_state);
              Cy_tmp2{i,j} = -eye(n_state);
          else
              Cy_tmp1{i,j} = zeros(n_state,n_state);
              Cy_tmp2{i,j} = zeros(n_state,n_state);
          end
       end
    end

    Cy = [cell2mat(Cy_tmp1);cell2mat(Cy_tmp2)];
end

author: ccalas

repo: mpc

save_path: github-repos/MATLAB/ccalas-mpc


text:
function [labelErrors, labelSharingErrors, wrongIndices] = findLabelErrors(realLabels, predLabels)

% FINDLABELERRORS Find the number of errors in the predicted vs the real
% labels
% COPYRIGHT: Andreas C. Damianou, 2012
% VARGPLVM

%%% WARNING: This function only works if there are no multi-label examples!

% First, transform any encoding to binary
% (e.g. if it is 1 -1 -1 -> 1 0 0)
lb = min(min(realLabels));
ub = max(max(realLabels));
realLabels(realLabels == lb) = 0;
realLabels(realLabels == ub) = 1;

predLabels(predLabels >= 0.5) = 1;
predLabels(predLabels < 0.5) = 0;

% %--- Old code
%     realUb = find(realLabels == ub);
%     predUb = predLabels(realUb);
%     % The number of times that a true label wasn't found
%     labelErrors = length(realUb) - sum(predUb); %equivalent to: length(find(~predUb))
%
 % The number of times more than one label was found

 labelSharingErrors = sum((sum(predLabels,2)>1));

 %--
 realLabels2 = transformLabels(realLabels);
 predLabels2 = transformLabels(predLabels);
 uniqueLabels = unique(realLabels2);
 labelErrors = 0;
 for j=1:length(uniqueLabels)
     trueClassInd = find(realLabels2 == uniqueLabels(j));
     wrongIndices{j} = find(~(predLabels2(trueClassInd) == uniqueLabels(j)));
     labelErrors = labelErrors + length(wrongIndices{j});
 end

author: SheffieldML

repo: vargplvm


text:
function [bhat sigmahatb sigmahat]=panel1estimates(X,Y,N,n,q,k,T)

% obtain an estimate of the VAR coefficients
% initiate the betahat vector
betahat=[];
% initiate the sum of the beta_i vectors and sigma_i matrices
betasum=zeros(q,1);
sigmasum=zeros(n,n);
% loop over units
for ii=1:N
% obtain Yi and Xi
Xi=X(:,:,ii);
Yi=Y(:,:,ii);
% estimate the VAR coefficients for this unit
Bhati=(Xi'*Xi)\(Xi'*Yi);
betahat(:,:,ii)=Bhati(:);
% estimate the residuals for this unit
EPShati=Yi-Xi*Bhati;
EPShat(:,:,ii)=EPShati;
% estimate the variance-covariance matrix for this unit and add to summation
sigmahati=(1/(T-k-1))*EPShati'*EPShati;
sigmasum=sigmasum+sigmahati;
% add the betahat value to the sum
betasum=betasum+betahat(:,:,ii);
end

% estimate bhat and sigmahat
bhat=(1/N)*betasum;
sigmahat=(1/N)*sigmasum;

% eventually estimate sigmab, the variance covariance matrix of the bhat vector of coefficients
sigmahatb=zeros(q,q);
for ii=1:N
sigmahatb=(betahat(:,1,ii)-bhat)*(betahat(:,1,ii)-bhat)';
end
sigmahatb=(1/(N*(N-1)))*sigmahatb;

author: european-central-bank

repo: BEAR-toolbox


text:
%
% Usage:   [val [paths]]=mexEvalPathCoding(U,DAG,param);
%
% Name: mexEvalPathCoding
%
% Description: mexEvalPathCoding evaluate the path coding penalies
%         of http://arxiv.org/abs/1204.4539 and provides a path
%         decomposition of a vector W.
%
%         Given an input matrix U=[u^1,\ldots,u^n],
%
%
% Inputs: U:  double p x n matrix   (input signals)
%               m is the signal size
%         DAG:  struct
%               with three fields, weights, start_weights, stop_weights
%         for a graph with |V| nodes and |E| arcs,
%         DAG.weights: sparse double |V| x |V| matrix. Adjacency
%               matrix. The non-zero entries represent costs on arcs
%               linking two nodes.
%         DAG.start_weights: dense double |V| vector. Represent the costs
%               of starting a path from a specific node.
%         DAG.stop_weights: dense double |V| vector. Represent the costs
%               of ending a path at a specific node.
%
%         if param.regul='graph-path-l0', non-convex penalty
%         if param.regul='graph-path-conv', convex penalty
%
%         param: struct
%               param.regul (choice of regularization, see above)
%               param.verbose (optional, verbosity level, false by default)
%               param.precision (optional, by default a very large integer.
%                 It returns approximate proximal operator by choosing a small integer,
%                 for example, 100 or 1000.
%               param.numThreads (optional, number of threads for exploiting
%                 multi-core / multi-cpus. By default, it takes the value -1,
%                 which automatically selects all the available CPUs/cores).
%
% Output: V: double 1 x n vector (values of the objective function)
%         paths: optional, double sparse p x k matrix. selected paths for the
%                first column of U
%
% Author: Julien Mairal, 2012

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function dz = rocketDynamics(z,u)
% dz = rocketDynamics(z,u)
%
% The basic dynamics and drag coefficient data are from the paper:
%
%   "Drag-law Effects in the Goddard Problem"
%   P. Tsiotras, H. Kelley, H.Kelley    1991
%
% INPUTS:
%   z = [3,n] = [h; v; m] = state vector
%   u = [1,n] = [T] = control = thrust
%

h = z(1,:);   %Height
v = z(2,:);   %Velocity
m = z(3,:);   %Mass
T = u;        %Thrust

%%%% Density of air:
% altitude = [0, 1e4, 2e4, 3e4, 4e4];  %(m)  %height above the ground
% density = [1.23, 0.41, 0.089, 0.018, 0.004];   %(kg/m^3)  density of air
density = 1.474085291*(0.9998541833.^h);  %Data fit off of wolfram alpha

%%%% Drag coefficient, calculated from paper:
A1 = 0.0095;
A2 = 25;
A3 = 0.953;
A4 = 0.036;
speedOfSound = 280;  %(m/s)  %At 10 km altitude
mach = v/speedOfSound;
Cd = A1*atan(A2*(mach-A3))+A4;

%%%% Compute the drag:
Area = pi*3.66;  %(m^2) cross-sectional area (SpaceX F9 Falcon)
D = 0.5*Cd.*Area.*density.*v.^2;

%%%% Compute gravity from inverse-square law:
rEarth = 6.3674447e6;  %(m) radius of earth
mEarth = 5.9721986e24;  %(kg) mass of earth
G = 6.67e-11; %(Nm^2/kg^2) gravitational constant
g = G*mEarth./((h+rEarth).^2);

%%%% Complete the calculation:
c = 4500;  %(m/s) rocket exhaust velocity - chemical rocket
dh = v;  %vertical velocity
dv = (T-D)./m - g;   %vertical acceleration
dm = -T/c;   %mass rate

dz = [dh;dv;dm];

end

author: MatthewPeterKelly

repo: OptimTraj


text:
function pbma_write_test ( file_out_name )

%*****************************************************************************80
%
%% PBMA_WRITE_TEST tests the ASCII portable bit map write routines.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 August 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character *file_out_name, the name of the file
%    to contain the ASCII portable pixel map data.
%
  nrow = 200;
  ncol = 200;
%
%  Set the data.
%
  b = pbma_example ( nrow, ncol );
%
%  Write the data to the file.
%
  pbma_write ( file_out_name, b );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function b = loadobj(a)
% loadobj for file_array class
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id$

if isa(a,'file_array')
    b = a;
else
    a = permission(a, 'rw');
    b = file_array(a);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function suite = test_assertElementsAlmostEqual
initTestSuite;

%===============================================================================
function test_happyCase

% All code here should execute with no error.
assertElementsAlmostEqual(1, 1 + sqrt(eps)/10);
assertElementsAlmostEqual(1, 1 + sqrt(eps)/10, 'custom message');

%===============================================================================
function test_failedAssertion

f = @() assertElementsAlmostEqual(1, 1 + 10*sqrt(eps));
assertExceptionThrown(f, 'assertElementsAlmostEqual:tolExceeded');

%===============================================================================
function test_nonFloatInputs()
assertExceptionThrown(@() assertElementsAlmostEqual('hello', 'world'), ...
    'assertElementsAlmostEqual:notFloat');

%===============================================================================
function test_sizeMismatch()
assertExceptionThrown(@() assertElementsAlmostEqual(1, [1 2]), ...
    'assertElementsAlmostEqual:sizeMismatch');

function test_finiteAndInfinite()
assertExceptionThrown(@() assertElementsAlmostEqual(1, Inf), ...
    'assertElementsAlmostEqual:tolExceeded');

function test_infiniteAndInfinite()
assertElementsAlmostEqual(Inf, Inf);

function test_finiteAndNaN()
assertExceptionThrown(@() assertElementsAlmostEqual(1, NaN), ...
    'assertElementsAlmostEqual:tolExceeded');

function test_nanAndNaN()
assertElementsAlmostEqual(NaN, NaN);

function test_plusMinusInfinity()
assertExceptionThrown(@() assertElementsAlmostEqual(+Inf, -Inf), ...
    'assertElementsAlmostEqual:tolExceeded');

function test_infiniteAndNaN()
assertExceptionThrown(@() assertElementsAlmostEqual(Inf, NaN), ...
    'assertElementsAlmostEqual:tolExceeded');

author: ricedsp


text:
function [L,a,b] = rgb2lab(R,G,B)
% function [L, a, b] = RGB2Lab(R, G, B)
% RGB2Lab takes matrices corresponding to Red, Green, and Blue, and
% transforms them into CIELab.  This transform is based on ITU-R
% Recommendation  BT.709 using the D65 white point reference.
% The error in transforming RGB -> Lab -> RGB is approximately
% 10^-5.  RGB values can be either between 0 and 1 or between 0 and 255.
% By Mark Ruzon from C code by Yossi Rubner, 23 September 1997.
% Updated for MATLAB 5 28 January 1998.

if (nargin == 1)
  B = double(R(:,:,3));
  G = double(R(:,:,2));
  R = double(R(:,:,1));
end

if ((max(max(R)) > 1.0) | (max(max(G)) > 1.0) | (max(max(B)) > 1.0))
  R = R/255;
  G = G/255;
  B = B/255;
end

[M, N] = size(R);
s = M*N;

% Set a threshold
T = 0.008856;

RGB = [reshape(R,1,s); reshape(G,1,s); reshape(B,1,s)];

% RGB to XYZ
MAT = [0.412453 0.357580 0.180423;
       0.212671 0.715160 0.072169;
       0.019334 0.119193 0.950227];
XYZ = MAT * RGB;

X = XYZ(1,:) / 0.950456;
Y = XYZ(2,:);
Z = XYZ(3,:) / 1.088754;

XT = X > T;
YT = Y > T;
ZT = Z > T;

fX = XT .* X.^(1/3) + (~XT) .* (7.787 .* X + 16/116);

% Compute L
Y3 = Y.^(1/3);
fY = YT .* Y3 + (~YT) .* (7.787 .* Y + 16/116);
L  = YT .* (116 * Y3 - 16.0) + (~YT) .* (903.3 * Y);

fZ = ZT .* Z.^(1/3) + (~ZT) .* (7.787 .* Z + 16/116);

% Compute a and b
a = 500 * (fX - fY);
b = 200 * (fY - fZ);

L = reshape(L, M, N);
a = reshape(a, M, N);
b = reshape(b, M, N);

if ((nargout == 1) | (nargout == 0))
  L = cat(3,L,a,b);
end

author: Cloud-CV

repo: object-proposals


text:
function vis_PlaneAxisProperties(ax,lattice,msize,pos)

% VIS_PLANEAXISPROPERTIES Set axis properties for SOM_CPLANE,
%                         SOM_PIEPLANE, SOM_BARPLANE and SOM_PLOTPLANE.
%
% vis_PlaneAxisProperties(ax,lattice,msize,pos)
%
%  Input arguments:
%   ax        (scalar) axis handle
%   lattice   (string) 'hexa', 'rect', 'hexaU' or 'rectU'
%             (matrix) defines the patch, see e.g. help vis_patch
%   msize     (vector) a 1x2 vector defining the grid size
%   pos       (vector) a 1x2 vector that determines position of
%                      origin or NaN which means default operation:
%                      origin to [1 1] and tighten axis limits
%                      according to the grid size.
%
% This is a subfunction for SOM_CPLANE, SOM_PIEPLANE, SOM_BARPLANE and
% SOM_PLOTPLANE. This subfunction sets the proper values for axis.

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

% Version 2.0beta Johan 060799

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xdim=msize(1);ydim=msize(2);
set(ax,'Visible','off');
set(get(ax,'Title'),'Visible','on');
set(ax,'XaxisLocation','Top');            % axis orientation
set(ax,'xdir','normal');                  % = axis ij = matrix mode
set(ax,'ydir','reverse');

switch lattice
case {'rect', 'rectU'}
  lelim=-.51; rilim=.51; uplim=-.51; lolim=.51;  % axis limits
  set(ax,'DataAspectRatio', [1 1 1]);            % =axis equal
case {'hexa','hexaU'}
  lelim=-.51; rilim=1.01; uplim=-.67; lolim=.67; % axis limits
  set(ax,'DataAspectRatio',[0.9015 1 1]);        % this corrects hexagons
end

% Nan: default origin [1 1] & tighten the axis
if isnan(pos)
  set(ax,'XLim',[1+lelim ydim+rilim],'YLim',[1+uplim xdim+lolim], ...
      'XLimMode','manual','YLimMode','manual'); % tighten the axis
end

author: ilarinieminen

repo: SOM-Toolbox


text:
function quatProd = quatMult(quat1, quat2)
%
% Multiplies two quaternions using the {i,j,k,1} convention
%
    if( size(quat1,1) ~= 4 || size(quat1,2) ~= 1 ...
        || size(quat2,1) ~= 4 || size(quat2,2) ~= 1 )
        error('Input quaternions must be 4x1');
    end

    quatProd = quatLeftComp(quat1) * quat2;
end

repo: MSCKF


text:
function [HDR] = leadidcodexyz(arg1)
% LeadIdCodeXYZ uses the Label information for computing the
%       LeadIdCode and the XYZ position of the EEG Electrodes
%       according to Annex A of FEF Vital Signs Format [1]
%
%   HDR = leadidcodexyz(HDR);
%	adds HDR.LeadIdCode and HDR.ELEC.XYZ, if needed.
%
% see also: SLOAD, SOPEN, PHYSICALUNITS, doc/leadidtable_scpecg.txt, doc/elecpos.txt
%
% Reference(s):
% [1] CEN/TC251/PT40 (2001)
% 	File Exchange Format for Vital Signs - Annex A
%
% Birbaumer, N. (2006). Brain-computer-interface research: Coming of age. Clinical Neurophysiology, 117:479–83.
%  http://www.acns.org/pdfs/ACFDD46.pdf.
% ACNS (2006). Guidelines for standard electrode position nomenclature. American Clinical
% Neurophysiology Society. http://www.acns.org/pdfs/ACFDD46.pdf.

%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 3
% of the License, or (at your option) any later version.

%	$Id$
%	Copyright (C) 2006,2007,2008,2009 by Alois Schloegl <a.schloegl@ieee.org>
%    	This is part of the BIOSIG-toolbox http://biosig.sf.net/

global BIOSIG_GLOBAL;
% BIOSIG_GLOBAL=[]; %%% used for debugging, only.

if ~isfield(BIOSIG_GLOBAL,'Phi');
	BIOSIG_GLOBAL.ISLOADED_XYZ = 0 ;
end;
if ~isfield(BIOSIG_GLOBAL,'ISLOADED_XYZ')
	BIOSIG_GLOBAL.ISLOADED_XYZ = 0 ;
end;
if ~BIOSIG_GLOBAL.ISLOADED_XYZ;
	f = which('getfiletype.m'); 	% identify path to biosig
	[p,f,e] = fileparts(f);
	[p,f,e] = fileparts(p);

        BIOSIG_GLOBAL.ISLOADED_XYZ = 0 ;

        N = 0;
        fid = fopen(fullfile(p,'doc','leadidtable_scpecg.txt'),'r');
        s = char(fread(fid,[1,inf],'uint8'));
        fclose(fid);

        Code = repmat(NaN, 200, 1); Phi = Code; Theta = Code;
        while ~isempty(s),
        	[t,s] = strtok(s,[10,13]);
                if ~length(t)
                elseif ~strncmp(t,'#',1)
                	ix3 = strfind(t,'MDC_ECG_LEAD_');
                	if isempty(ix3)
                		ix3 = length(t)+1;
                	end
                        [t1,t2] = strtok(t(1:ix3-1),[9,32]);
                        [t2,t3] = strtok(t2,[9,32]);
                        id = str2double(t2);
                        N  = N + 1;
                        Labels{N,1}	  = t1;
                        Code(N,1)         = id;
                        Description{N,1}  = deblank(t3);
                        %MDC_ECG_LEAD{N,1} = t(ix3+13:end)
                        MDC_ECG_LEAD{N,1} = t(ix3:end);
                end;
        end;
        N1 = N;

        % load table
        fid = fopen(fullfile(p,'doc','elecpos.txt'),'r');
        t = char(fread(fid,[1,inf],'uint8'));
        fclose(fid);

        % extract table information
        while ~isempty(t)
                [x,t] = strtok(t,[10,13]);
                if isempty(x)
                elseif strncmp(x,'#',1)
                else
                        N = N + 1;
                        [num,status,strarray] = str2double(x);
                        Code(N,1)   = num(1);
                        Labels{N,1} = upper(strarray{2});
                        Phi(N,1)    = num(3);
                        Theta(N,1)  = num(4);
                end;
        end;
        Phi   = Phi(:)  *pi/180;
        Theta = Theta(:)*pi/180;

        % loading is done only once.
        BIOSIG_GLOBAL.XYZ = [sin(Theta).*cos(Phi), sin(Theta).*sin(Phi), cos(Theta)];
        BIOSIG_GLOBAL.Phi          = Phi*180/pi;
        BIOSIG_GLOBAL.Theta        = Theta*180/pi;
        BIOSIG_GLOBAL.LeadIdCode   = Code;
        BIOSIG_GLOBAL.Label        = Labels;
        BIOSIG_GLOBAL.Description  = Description;
        BIOSIG_GLOBAL.MDC_ECG_LEAD = MDC_ECG_LEAD;

        BIOSIG_GLOBAL.ISLOADED_XYZ = 1;
end;

if nargin<1,
        HDR.LeadIdCode = BIOSIG_GLOBAL.LeadIdCode;
        HDR.Label      = BIOSIG_GLOBAL.Label;
        HDR.ELEC.XYZ   = BIOSIG_GLOBAL.XYZ;
        HDR.TYPE       = 'ELPOS';

else    % electrode code and position

        if isstruct(arg1)
                HDR = arg1;
        elseif isnumeric(arg1),
                HDR.LeadIdCode = arg1;
        else
                HDR.Label = arg1;
        end;

        tmp.flag1 = isfield(HDR,'ELEC');
        if tmp.flag1,
                tmp.flag1 = isfield(HDR.ELEC,'XYZ');
        end;
        if tmp.flag1,
                tmp.flag1 = any(HDR.ELEC.XYZ(:));
        end;
        tmp.flag2 = isfield(HDR,'LeadIdCode');
        tmp.flag3 = isfield(HDR,'Label');

        if (~tmp.flag1 || ~tmp.flag2 || ~tmp.flag3),
        	if 0,
        	elseif tmp.flag3,
                        if ischar(HDR.Label)
                                HDR.Label = cellstr(HDR.Label);
                        end;
	                NS = length(HDR.Label);
	        elseif tmp.flag2,
	        	NS = length(HDR.LeadIdCode);
        	elseif isfield(HDR,'NS')
        		NS = HDR.NS;
        		HDR.LeadIdCode = zeros(1,HDR.NS);
	        end;

                if tmp.flag3,
	                if ~tmp.flag1,
				HDR.ELEC.XYZ   = repmat(NaN,NS,3);
				HDR.ELEC.Phi   = repmat(NaN,NS,1);
				HDR.ELEC.Theta = repmat(NaN,NS,1);
	        	end;
                	if ~tmp.flag2,
                       		HDR.LeadIdCode = repmat(NaN,NS,1);
	        	end;
                	for k = 1:NS;
				Label = upper(deblank(HDR.Label{k}));
				pos = find(Label==':');
				if ~isempty(pos)
					Label = Label(pos+1:end);
                            	end;
	                        ix = strmatch(Label,BIOSIG_GLOBAL.Label,'exact');

	                        if length(ix)==2,
	                        	%%%%% THIS IS A HACK %%%%%
	                        	%% solve ambiguity for 'A1','A2'; could be EEG or ECG
	                        	if sum(HDR.LeadIdCode(1:k)>=996)>sum(HDR.LeadIdCode(1:k)<996)
	                        		%% majority are EEG electrodes,
	                        		ix = ix(find(BIOSIG_GLOBAL.LeadIdCode(ix)>996));
	                        	else
	                        		%% majority are ECG electrodes,
	                        		ix = ix(find(BIOSIG_GLOBAL.LeadIdCode(ix)<996));
	                        	end;
	                        elseif isempty(ix)
		                        ix = strmatch(deblank(HDR.Label{k}),BIOSIG_GLOBAL.MDC_ECG_LEAD,'exact');
	                        end;

        	                if (length(ix)==1),
	                	        if ~tmp.flag1,
        	                	        HDR.ELEC.XYZ(k,1:3) = BIOSIG_GLOBAL.XYZ(ix,:);
                	        	        HDR.ELEC.Phi(k)   = BIOSIG_GLOBAL.Phi(ix);
                        		        HDR.ELEC.Theta(k) = BIOSIG_GLOBAL.Theta(ix);
                        		end;
                	        	if ~tmp.flag2,
                                		HDR.LeadIdCode(k,1) = BIOSIG_GLOBAL.LeadIdCode(ix);
	                        	end;
        	                end;
                        end;
                else
			HDR.Label = cell(NS,1);
			for k = 1:NS;
				ix = find(BIOSIG_GLOBAL.LeadIdCode==HDR.LeadIdCode(k));
				if (length(ix)>=1),
					ix=ix(1);	%
	                                HDR.Label{k} = BIOSIG_GLOBAL.Label{ix};
					if ~tmp.flag1,
        		                        HDR.ELEC.XYZ(k,1:3) = BIOSIG_GLOBAL.XYZ(ix,1:3);
        		                        HDR.ELEC.Phi(k,1)   = BIOSIG_GLOBAL.Phi(ix);
        		                        HDR.ELEC.Theta(k,1) = BIOSIG_GLOBAL.Theta(ix);
					end;
				else
					HDR.Label{k} = ['#',int2str(k)];
				end;
			end;
                end;
	end;
	if tmp.flag3 && ~any(HDR.LeadIdCode),
               	for k = 1:HDR.NS;
			ix = strmatch(upper(HDR.Label{k}),BIOSIG_GLOBAL.Label,'exact');
			if length(ix)==1,
				HDR.LeadIdCode(k) = BIOSIG_GLOBAL.LeadIdCode(ix);
			end;
		end;
        end;
end;


text:
function [alpha,beta]=cnonlin2(evalstr,x,y,sy,p,delta,lambda);

% CNONLIN2 is a subsidiary function of CNONLIN. It must be called through CNONLIN.

m=length(p);
beta=-0.5*firstd(evalstr,x,y,sy,p,delta);
d2XiSq_da2=zeros(m);
for j=1:m;
  p(j)=p(j)+0.5*delta(j);
  dXiSqPlus=firstd(evalstr,x,y,sy,p,delta/1000);
  p(j)=p(j)-delta(j);
  dXiSqMinus=firstd(evalstr,x,y,sy,p,delta/1000);
  p(j)=p(j)+0.5*delta(j);
  d2XiSq_da2(:,j)=(dXiSqPlus-dXiSqMinus)/delta(j);
 end;
alpha=0.5*d2XiSq_da2;
for i=1:m
 alpha(i,i)=alpha(i,i)*(1+lambda);
end;

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function  SELVEparam = initSELVE(traindata, SELVEparam)

SELVEparam.s = 2;   %number of nearest anchors, please tune this parameter on different datasets
SELVEparam.lambda = 0.1;
SELVEparam.beta = 0.5;
SELVEparam.RedDim = 300;  % the left dimensions of original data
SELVEparam.m = 300;       % the number of landmark
SELVEparam.sigma = 0;

kmMaxIter = 10;
kmNumRep = 1;
[label,anchor] = litekmeans(traindata,SELVEparam.m,'MaxIter',kmMaxIter,'Replicates',kmNumRep);
clear kmMaxIter kmNumRep;

temptraindata = zeros(size(traindata,1),size(traindata,2));
tempanchor = zeros(SELVEparam.m,size(traindata,2));
totalnum = 0;

for i=1:length(unique(label))
    tempnum = find(label == i);
    temptraindata(totalnum +1: totalnum + length(tempnum),:) = traindata(tempnum,:);
    tempanchor(totalnum +1:totalnum + length(tempnum),:) = repmat(anchor(i,:),length(tempnum),1);
    totalnum = totalnum + length(tempnum);
end
Gamma = temptraindata - tempanchor;  % ins * fea
SELVEparam.GtG = Gamma'*Gamma;
clear temptraindata tempanchor Gamma tempnum totalnum;

%% searching for transformation matrix
SELVEparam.XtX = traindata'*traindata;
SELVEparam.anchor = anchor;
SELVEparam.label = label;

author: willard-yuan

repo: hashing-baseline-for-image-retrieval


text:
function pde = elli3DlinIntf2(bm,bp,delta)
%% USAGE: polynomial solution for Poisson equation
%  Last Modified: 12/15/2020 by GRC

%% PDE Structure
pde = struct('intf',@intf,'f',@f,'fm',@fm,'fp',@fp,'exactu',@exactu,...
    'um',@um,'up',@up,'Dxu',@Dxu,'Dxum',@Dxum,'Dxup',@Dxup,'Dyu',@Dyu,...
    'Dyum',@Dyum,'Dyup',@Dyup,'Dzu',@Dzu,'Dzum',@Dzum,'Dzup',@Dzup,...
    'A',@A,'Am',@Am,'Ap',@Ap,'one',@one,'gD',@gD);

pde.bm = bm;
pde.bp = bp;
%% interface function
    function u = intf(x,y,z)
        u = (x-delta);
    end

%% exact solution
    function u = exactu(x,y,z)
        u = um(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up(x(id),y(id),z(id));
    end
    function u = um(x,y,z)
        r = (1+delta)*sin((x-delta)*pi/(1+delta)).*sin(pi*y).*sin(pi*z);
        u = r/bm;
    end
    function u = up(x,y,z)
        r = (1-delta)*sin((x-delta)*pi/(1-delta)).*sin(pi*y).*sin(pi*z);
        u = r/bp;
    end
%% Boundary Function
    function u = gD(x,y,z)
        u = exactu(x,y,z);
    end
%% Derivative of the exact solution
    function u = Dxu(x,y,z)
        u = Dxum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dxup(x(id),y(id),z(id));
    end
    function u = Dxum(x,y,z)
        r = cos((x-delta)*pi/(1+delta)).*sin(pi*y).*sin(pi*z)*pi;
        u = r/bm;
    end
    function u = Dxup(x,y,z)
        r = cos((x-delta)*pi/(1-delta)).*sin(pi*y).*sin(pi*z)*pi;
        u = r/bp;
    end

    function u = Dyu(x,y,z)
        u = Dyum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dyup(x(id),y(id),z(id));
    end
    function u = Dyum(x,y,z)
        r = (1+delta)*sin((x-delta)*pi/(1+delta)).*cos(pi*y).*sin(pi*z)*pi;
        u = r/bm;
    end
    function u = Dyup(x,y,z)
        r = (1-delta)*sin((x-delta)*pi/(1-delta)).*cos(pi*y).*sin(pi*z)*pi;
        u = r/bp;
    end

    function u = Dzu(x,y,z)
        u = Dzum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dzup(x(id),y(id),z(id));
    end
    function u = Dzum(x,y,z)
        r = (1+delta)*sin((x-delta)*pi/(1+delta)).*sin(pi*y).*cos(pi*z)*pi;
        u = r/bm;
    end
    function u = Dzup(x,y,z)
        r = (1-delta)*sin((x-delta)*pi/(1-delta)).*sin(pi*y).*cos(pi*z)*pi;
        u = r/bp;
    end

%% right hand side function
    function u = f(x,y,z)
        u = fm(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp(x(id),y(id),z(id));
    end
    function u = fm(x,y,z)
        r = sin((x-delta)*pi/(1+delta)).*sin(pi*y).*sin(pi*z)*pi^2;
        u = r*(1/(1+delta)+2*(1+delta));
    end
    function u = fp(x,y,z)
        r = sin((x-delta)*pi/(1-delta)).*sin(pi*y).*sin(pi*z)*pi^2;
        u = r*(1/(1-delta)+2*(1-delta));
    end

%% Diffusion coefficient function
    function u = A(x,y,z)
        u = Am(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Ap(x(id),y(id),z(id));
    end
    function u = Am(x,y,z)
        u = bm*ones(size(x));
    end
    function u = Ap(x,y,z)
        u = bp*ones(size(x));
    end

%% Other function
    function u = one(x,y,z)
        u = ones(size(x));
    end
end

repo: ifem


text:
function showrateh4(h1,err1,k1,opt1,str1,h2,err2,k2,opt2,str2,...
                    h3,err3,k3,opt3,str3,h4,err4,k4,opt4,str4)
%% SHOWRATEH4 rate of two error sequences
%
% showrateh4(N1,err1,k1,opt1,str1,N2,err2,k2,opt2,str2,N3,err3,k3,opt3,str3,err4,k4,opt4,str4)
% plots the err1 vs N1, err2 vs N2, err3 vs N3, and err4 vs N4 in the
% loglog scale. Additional input
%
%   - k1, k2, k3, k4: specify the starting indices; see showrate
%   - opt1, opt2, opt3, opt4: the line color and style
%   - str1, str2, str3, str4: strings used in legend
%
% Example
%
%
% See also showrate, showresult, showmesh, showsolution
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

N1 = 1./h1; N2 = 1./h2; N3 =1./h3; N4 =1./h4;
if (nargin<=2)
    k1 = 1; opt1 = '-*';
end
r1 = showrate(N1,err1,k1,opt1);
hold on
r2 = showrate(N2,err2,k2,opt2);
r3 = showrate(N3,err3,k3,opt3);
r4 = showrate(N4,err4,k4,opt4);
h_legend = legend(str1,['C_1h^{' num2str(-r1) '}'],...
                  str2,['C_2h^{' num2str(-r2) '}'],...
                  str3,['C_3h^{' num2str(-r3) '}'],...
                  str4,['C_4h^{' num2str(-r4) '}'],'LOCATION','Best');
set(h_legend,'FontSize',12);
xlabel('log(1/h)');

repo: Aero-matlab


text:
function [ row, col, a, seed ] = wathen_st ( nx, ny, nz_num, seed )

%*****************************************************************************80
%
%% WATHEN_ST: Wathen matrix stored in sparse triplet format.
%
%  Discussion:
%
%    When dealing with sparse matrices in MATLAB, it can be much more efficient
%    to work first with a triple of I, J, and X vectors, and only once
%    they are complete, convert to MATLAB's sparse format.
%
%    The Wathen matrix is a finite element matrix which is sparse.
%
%    The entries of the matrix depend in part on a physical quantity
%    related to density.  That density is here assigned random values between
%    0 and 100.
%
%    The matrix order N is determined by the input quantities NX and NY,
%    which would usually be the number of elements in the X and Y directions.
%
%    The value of N is
%
%      N = 3*NX*NY + 2*NX + 2*NY + 1,
%
%    The matrix is the consistent mass matrix for a regular NX by NY grid
%    of 8 node serendipity elements.
%
%    The local element numbering is
%
%      3--2--1
%      |     |
%      4     8
%      |     |
%      5--6--7
%
%    Here is an illustration for NX = 3, NY = 2:
%
%     23-24-25-26-27-28-29
%      |     |     |     |
%     19    20    21    22
%      |     |     |     |
%     12-13-14-15-16-17-18
%      |     |     |     |
%      8     9    10    11
%      |     |     |     |
%      1--2--3--4--5--6--7
%
%    For this example, the total number of nodes is, as expected,
%
%      N = 3 * 3 * 2 + 2 * 2 + 2 * 3 + 1 = 29
%
%    The matrix is symmetric positive definite for any positive values of the
%    density RHO(X,Y).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 July 2014
%
%  Author:
%
%    Original MATLAB version by Nicholas Higham.
%    Modifications by Tim Davis.
%    Modifications by John Burkardt.
%
%  Reference:
%
%    Nicholas Higham,
%    Algorithm 694: A Collection of Test Matrices in MATLAB,
%    ACM Transactions on Mathematical Software,
%    Volume 17, Number 3, September 1991, pages 289-305.
%
%    Andrew Wathen,
%    Realistic eigenvalue bounds for the Galerkin mass matrix,
%    IMA Journal of Numerical Analysis,
%    Volume 7, Number 4, October 1987, pages 449-457.
%
%  Parameters:
%
%    Input, integer NX, NY, values which determine the size of the matrix.
%
%    Input, integer NZ_NUM, the number of values used to describe the matrix.
%
%    Input/output, integer SEED, the random number seed.
%
%    Output, integer ROW(NZ_NUM), COL(NZ_NUM), the row and column indices
%    of the nonzero entries.
%
%    Output, real A(NZ_NUM), the nonzero values.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'WATHEN_ST - Fatal error!\n' );
    fprintf ( 1, '  Not enough input.\n' );
    error ( 'WATHEN_ST - Fatal error!' );
  end

  if ( nargin < 2 )
    ny = nx;
    fprintf ( 1, '\n' );
    fprintf ( 1, '  NY was not supplied.  Setting NY = NX = %d.\n', ny );
  end

  if ( nargin < 3 )
    nz_num = wathen_st_size ( nx, ny );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  NZ_NUM was not supplied.  NZ_NUM = %d\n', nz_num );
  end

  if ( nargin < 4 )
    seed = 123456789;
  end

  em = [
     6.0, -6.0,  2.0, -8.0,  3.0, -8.0,  2.0, -6.0;
    -6.0, 32.0, -6.0, 20.0, -8.0, 16.0, -8.0, 20.0;
     2.0, -6.0,  6.0, -6.0,  2.0, -8.0,  3.0, -8.0;
    -8.0, 20.0, -6.0, 32.0, -6.0, 20.0, -8.0, 16.0;
     3.0, -8.0,  2.0, -6.0,  6.0, -6.0,  2.0, -8.0;
    -8.0, 16.0, -8.0, 20.0, -6.0, 32.0, -6.0, 20.0;
     2.0, -8.0,  3.0, -8.0,  2.0, -6.0,  6.0, -6.0;
    -6.0, 20.0, -8.0, 16.0, -8.0, 20.0, -6.0, 32.0 ]';

  row = zeros(nz_num,1);
  col = zeros(nz_num,1);
  a = zeros(nz_num,1);

  node = zeros(8,1);

  k = 0;

  for j = 1 : ny
    for i = 1 : nx

      node(1) = 3 * j * nx + 2 * i + 2 * j + 1;
      node(2) = node(1) - 1;
      node(3) = node(2) - 1;
      node(4) = ( 3 * j - 1 ) * nx + 2 * j + i - 1;
      node(5) = 3 * ( j - 1 ) * nx + 2 * i + 2 * j - 3;
      node(6) = node(5) + 1;
      node(7) = node(6) + 1;
      node(8) = node(4) + 1;

      [ rho, seed ] = r8_uniform_01 ( seed );
      rho = 100.0 * rho;

      for krow = 1 : 8
        for kcol = 1 : 8
          k = k + 1;
          row(k) = node(krow);
          col(k) = node(kcol);
          a(k) = rho * em(krow,kcol);
        end
      end

    end
  end

  return
end


text:
function I = gif2RGB(gifFile)
[x, map] = imread(gifFile);
I = ind2rgb(x, map);

author: Sable

repo: mcbench-benchmarks


text:
% set default value for a field if not set yet.
function value = ReturnFieldWithDefaultValue(para, fieldName, defaultValue)
    if ~isfield(para, fieldName)
        value = defaultValue;
    else
        value = para.(fieldName);
    end
end

author: singaxiong

repo: SignalGraph


text:
% VL_DEMO_DSIFT Demo: DSIFT

randn('state',0) ;
rand('state',0) ;

% read a test image
I = imread(fullfile(vl_root,'data','roofs1.jpg')) ;
I = single(vl_imdown(rgb2gray(I))) ;

% --------------------------------------------------------------------
%                                                      Basic benchmark
% --------------------------------------------------------------------

binSize = 4 ; % bin size in pixels
magnif = 3 ; % bin size / keypoint scale

elaps_dsift = [] ;
elaps_dsift_fast = [] ;
err_dsift = [] ;
err_dsift_fast = [] ;

binSizeRange = [3 4 5 6] ;
for wi = 1:length(binSizeRange)
  binSize = binSizeRange(wi) ;
  scale = binSize / magnif ;

  tic ;
  [f, d] = vl_dsift(vl_imsmooth(I, sqrt(scale.^2 - .25)), ...
                    'size', binSize, ...
                    'step', 2, ...
                    'bounds', [20,20,210,140], ...
                    'floatdescriptors', ...
                    'verbose') ;
  elaps_dsift(wi) = toc ;

  tic ;
  [f, dfast] = vl_dsift(vl_imsmooth(I, sqrt(scale.^2 - .25)), ...
                        'size', binSize, ...
                        'step', 2, ...
                        'bounds', [20,20,210,140], ...
                        'floatdescriptors', ...
                        'fast', ...
                        'verbose') ;
  elaps_dsift_fast(wi) = toc ;

  numKeys = size(f, 2) ;
  f_ = [f ; ones(1, numKeys) * scale ; zeros(1, numKeys)] ;

  tic ;
  [f_, d_] = vl_sift(I, ...
                     'magnif', magnif, ...
                     'frames', f_, ...
                     'firstoctave', -1, ...
                     'levels', 5, ...
                     'floatdescriptors') ;
  elaps_sift(wi) = toc ;

  err_dsift(wi)      = mean(mean(abs(d     - d_)) ./ mean(d_)) * 100 ;
  err_dsift_fast(wi) = mean(mean(abs(dfast - d_)) ./ mean(d_)) * 100 ;
end

figure(1) ; clf ; title('Descriptor  SIFT') ;
plot(binSizeRange, [err_dsift ; err_dsift_fast]', 'linewidth', 3) ;
legend('DSIFT', 'DSIFT fast') ;
ylabel('Approx error (%)') ;
xlabel('binSize parameter') ;
grid on ;
axis square ;

figure(2) ; title('Speedup on regular SIFT') ;
plot(binSizeRange, ...
     [elaps_sift ./ elaps_dsift ; ...
      elaps_sift ./ elaps_dsift_fast ; ...
      elaps_sift ./ elaps_sift], 'linewidth', 3) ;
legend('DSIFT', 'DSIFT fast', 'SIFT') ;
ylabel('Speedup') ;
xlabel('binSize parameter') ;
grid on ;
axis square ;

figure(1) ; vl_demo_print('dsift_accuracy') ;
figure(2) ; vl_demo_print('dsift_speedup') ;


text:
function [val1,val2,val3,val4] = outtowks(mode)
%
% Utility Function: OUTTOWKS
%
% This function provides a gateway from the environment to the workspace

% Author: Craig Borghesani
% Date: 9/5/94
% Revised: 11/15/94
% Copyright (c) 1999, Prentice-Hall

if mode == 0,
 f = gcf;
else
 f2 = gcf;
 f = get(f2,'userdata');
end
ui_data = get(f,'userdata');
ui_han = ui_data{1};
ui_obj = ui_data{2};
ui_han2 = get(ui_han(51),'userdata');
cur_sys = get(ui_han(30),'userdata');
stat_bar = get(ui_han(43),'userdata');
win_str = get(ui_han(48),'userdata');
if length(ui_han2),
 stat_bar2 = ui_han2(16);
end

if mode==0, % Setup

 if ~figflag(['Out to Workspace...',win_str]),
  f2 = figure('numbertitle','off','menubar','none','name',['Out to Workspace...',win_str],...
              'pos',[30,50,395,187],'color',[0.8,0.8,0.8],'resize','off',...
              'userdata',f);

% Variable Input Fields
  del=17;
  uicontrol('style','text','string','Variable Names',...
            'pos',[10,145+del,375,17],'horizontalalignment','left',...
            'backgroundcolor',[0.8,0.8,0.8]);
  ui_han2(1) = uicontrol('style','text','string','Num:',...
            'pos',[0,123+del,40,17],'horizontalalignment','right',...
            'backgroundcolor',[0.8,0.8,0.8]);
  ui_han2(2) = uicontrol('style','edit','pos',[40,123+del,55,20],'backgroundcolor','w',...
                     'callback','outtowks(1)',...
                     'horiz','left');

  ui_han2(3) = uicontrol('style','text','string','Den:',...
            'pos',[95,123+del,40,17],'horizontalalignment','right',...
            'backgroundcolor',[0.8,0.8,0.8]);
  ui_han2(4) = uicontrol('style','edit','pos',[135,123+del,55,20],'backgroundcolor','w',...
                     'callback','outtowks(1)',...
                     'horiz','left');

  ui_han2(5) = uicontrol('style','text','string','',...
            'pos',[190,123+del,40,17],'horizontalalignment','right','vis','off',...
            'backgroundcolor',[0.8,0.8,0.8]);
  ui_han2(6) = uicontrol('style','edit','pos',[230,123+del,55,20],'backgroundcolor','w',...
                     'callback','outtowks(1)','vis','off',...
                     'horiz','left');

  ui_han2(7) = uicontrol('style','text','string','',...
            'pos',[285,123+del,40,17],'horizontalalignment','right','vis','off',...
            'backgroundcolor',[0.8,0.8,0.8]);
  ui_han2(8) = uicontrol('style','edit','pos',[325,123+del,55,20],'backgroundcolor','w',...
                     'callback','outtowks(1)','vis','off',...
                     'horiz','left');

% Option Fields
  uicontrol('style','text','string','Format',...
            'pos',[10,100+del,180,17],'horizontalalignment','left',...
            'backgroundcolor',[0.8,0.8,0.8]);
  ui_han2(9) = uicontrol('style','radio','pos',[10,80+del,180,20],...
                     'string','Numerator/Denominator',...
                     'backgroundcolor',[0.8,0.8,0.8],...
                     'callback','outtowks(2)','value',1);
  ui_han2(10) = uicontrol('style','radio','pos',[10,60+del,180,20],...
                     'string','Zero/Pole/Gain',...
                     'backgroundcolor',[0.8,0.8,0.8],...
                     'callback','outtowks(2)','value',0);
  ui_han2(11) = uicontrol('style','radio','pos',[10,40+del,180,20],...
                     'string','State Space',...
                     'backgroundcolor',[0.8,0.8,0.8],...
                     'callback','outtowks(2)','value',0);

  uicontrol('style','text','string','Transfer Function',...
            'pos',[205,100+del,180,17],'horizontalalignment','left',...
            'backgroundcolor',[0.8,0.8,0.8]);
  ui_han2(12) = uicontrol('style','radio','pos',[205,80+del,180,20],...
                     'string','Plant',...
                     'backgroundcolor',[0.8,0.8,0.8],...
                     'callback','outtowks(3)','value',1);
  ui_han2(13) = uicontrol('style','radio','pos',[205,60+del,180,20],...
                     'string','Forward Controller',...
                     'backgroundcolor',[0.8,0.8,0.8],...
                     'callback','outtowks(3)','value',0);
  ui_han2(14) = uicontrol('style','radio','pos',[205,40+del,180,20],...
                     'string','Feedback Controller',...
                     'backgroundcolor',[0.8,0.8,0.8],...
                     'callback','outtowks(3)','value',0);

  for c=1:3,
   set(ui_han2(8+c),'userdata',ui_han2(:,8+[1:(c-1),(c+1):3]));
  end

  for c=1:3,
   set(ui_han2(11+c),'userdata',ui_han2(:,11+[1:(c-1),(c+1):3]));
  end

  ui_han2(15) = uicontrol('style','push','pos',[103,10+del,60,20],'string','OK');
  uicontrol('style','push','pos',[166,10+del,60,20],'string','Cancel',...
            'callback','outtowks(5)');
  uicontrol('style','push','pos',[232,10+del,60,20],'string','Help',...
            'callback','userhelp(5,3)');

  ui_han2(16) = uicontrol('style','edit','pos',[0,0,395,20],...
             'foregroundcolor','k',...
             'horizontalalignment','left');
  stat_bar2 = ui_han2(16);
  set(ui_han(51),'userdata',ui_han2);
 end

 set(stat_bar2,'string','Enter variable names into edit boxes, press OK when done');
 drawnow;

elseif mode==1, % edit fields

% build output and input strings
 output_str = [];

 for k = 2:2:8,
  var_str = get(ui_han2(k),'string');
  if length(var_str),
   output_str = [output_str,var_str,','];
  end
 end

% remove extra comma
 output_str(length(output_str))=[];

% setup callback to pass desired variables to workspace
 callback_str = ['[',output_str,']=outtowks(4);'];
 set(ui_han2(15),'callback',callback_str);

elseif mode==2, % format

 cur_obj = get(f2,'currentobject');
 data = get(cur_obj,'userdata');
 set(data,'value',0);
 set(cur_obj,'value',1);

 if cur_obj == ui_han2(9), % num/den
  set(ui_han2(5:8),'vis','off');
  set(ui_han2(1),'string','Num:');
  set(ui_han2(3),'string','Den:');

 elseif cur_obj == ui_han2(10), % z/p/k
  set(ui_han2(7:8),'vis','off');
  set(ui_han2(5:6),'vis','on');
  set(ui_han2(1),'string','Zero:');
  set(ui_han2(3),'string','Pole:');
  set(ui_han2(5),'string','Gain:');

 elseif cur_obj == ui_han2(11), % state space
  set(ui_han2(5:8),'vis','on');
  set(ui_han2(1),'string','A:');
  set(ui_han2(3),'string','B:');
  set(ui_han2(5),'string','C:');
  set(ui_han2(7),'string','D:');
 end

elseif mode==3, % plant or controller

 cur_obj = get(f2,'currentobject');
 data = get(cur_obj,'userdata');
 set(data,'value',0);
 set(cur_obj,'value',1);

elseif mode==4, % OK

 if get(ui_han2(12),'value'), sel = 3;
 elseif get(ui_han2(13),'value'), sel = 6;
 elseif get(ui_han2(14),'value'), sel = 9; end
 term_mat = get(ui_han(sel+cur_sys),'userdata');

 if get(ui_han2(9),'value'),
  [var1,var2]=termextr(term_mat);
  len = 2;
 end

 if get(ui_han2(10),'value'),
  [var1,var2,var3]=term2zpk(term_mat);
  len = 3;
 end

 if get(ui_han2(11),'value'),
  [n,d]=termextr(term_mat);
  [var1,var2,var3,var4]=tf2ss(n,d);
  len = 4;
 end

 ct=1;
 for k = 1:len,
  eval(['val',int2str(ct),'=var',int2str(k),';']);
  ct=ct+1;
 end

 set(stat_bar,'string','Variables passed to workspace');

 set(f2,'vis','off');

elseif mode==5, % Cancel

 set(f2,'vis','off');

end


text:
function [sys,x0,str,ts] = SimpleVehicleSFunction(t,x,u,flag)
% This file is a s-function template for simulating the simple vehicle
% model in Simulink.

% Choosing tire model
TireModel = VehicleDynamicsLateral.TirePacejka();
% Defining tire parameters
TireModel.a0        = 1;
TireModel.a1        = 0;
TireModel.a2        = 800;
TireModel.a3        = 3000;
TireModel.a4        = 50;
TireModel.a5        = 0;
TireModel.a6        = 0;
TireModel.a7        = -1;
TireModel.a8        = 0;
TireModel.a9        = 0;
TireModel.a10       = 0;
TireModel.a11       = 0;
TireModel.a12       = 0;
TireModel.a13       = 0;

% Choosing vehicle model
VehicleModel = VehicleDynamicsLateral.VehicleSimpleNonlinear();
% Defining vehicle parameters
VehicleModel.mF0    = 700;
VehicleModel.mR0    = 600;
VehicleModel.IT     = 10000;
VehicleModel.lT     = 3.5;
VehicleModel.nF     = 2;
VehicleModel.nR     = 2;
VehicleModel.wT     = 2;
VehicleModel.muy    = 0.8;
VehicleModel.tire   = TireModel;

switch flag

  %%%%%%%%%%%%%%%%%%
  % Initialization %
  %%%%%%%%%%%%%%%%%%
  case 0
    [sys,x0,str,ts]=mdlInitializeSizes();

  %%%%%%%%%%%%%%%
  % Derivatives %
  %%%%%%%%%%%%%%%
  case 1
    sys=mdlDerivatives(t,x,u,VehicleModel);

  %%%%%%%%%%%
  % Outputs %
  %%%%%%%%%%%
  case 3
    sys=mdlOutputs(t,x,u,VehicleModel);

  %%%%%%%%%%%%%%%%%%%
  % Unhandled flags %
  %%%%%%%%%%%%%%%%%%%
  case { 2, 4, 9 }
    sys = [];

  %%%%%%%%%%%%%%%%%
  % Vehicle model %
  %%%%%%%%%%%%%%%%%
  % Case 5 returns the vehicle model.
  case 5
    sys = VehicleModel;
    x0  =1; % Dummy
    str =1; % Dummy
    ts  =1; % Dummy

  %%%%%%%%%%%%%%%%%%%%
  % Unexpected flags %
  %%%%%%%%%%%%%%%%%%%%
  otherwise
    DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));

end
% end csfunc

%
%=============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%

function [sys,x0,str,ts]=mdlInitializeSizes()

% Definitions
sizes = simsizes;
sizes.NumContStates  = 6;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 6;
sizes.NumInputs      = 4;
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

% Setting initial conditions
VEL0 = 50/3.6; % Initial velocity [m/s]

x0  = [0 0 0 VEL0 0 0];
str = [];
ts  = [0 0];

% end mdlInitializeSizes
%
%=============================================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%=============================================================================
%

function sys = mdlDerivatives(t,x,u,vehicle)

% Defining input
vehicle.deltaf  = u(1);
vehicle.deltar  = u(2);
vehicle.Fxf     = u(3);
vehicle.Fxr     = u(4);

% Getting the vehicle model function (state equations)
ModelFunction   = @vehicle.Model;

sys = ModelFunction(t,x,0);

% end mdlDerivatives
%
%=============================================================================
% mdlOutputs
% Return the block outputs.
%=============================================================================
%

function sys=mdlOutputs(~,x,~,~)

% Output are all state variables
sys = x;

% end mdlOutputs

author: andresmendes

repo: Vehicle-Dynamics-Lateral


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: density estimation using a histogram
%
%==============================================================================

clear, close all, help(mfilename);

n = 1000;
x = linspace(0,2*pi,n); T = 0.5*(sin(x)+1); % discretized function
minT = 0; maxT = 1;                         % bounds for the bins
numberBins = 5;                             % number of bins
binWidth   = (maxT-minT)/numberBins;        % bin width
bins       = 0:binWidth:maxT;               % the bins
binsExt    = [-inf,bins(2:end-1),inf];      % don't miss anything

rhoHat = histc(T,binsExt); % compute the histogram and plot it
bar(bins+binWidth/2,rhoHat,0.99,'edgecolor','w','facecolor',0.8*[1,1,1]);
axis([minT,maxT,0,inf])

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


text:
function [p q m1 n1 tol] = spqr_singletons (A,tol)                          %#ok
%SPQR_SINGLETONS finds the singleton permutation of a sparse matrix A.
% [p q m1 n1 tol] = spqr_singletons (A, tol) finds permutation vectors p and q
% so that C=A(p,q) is a block 2-by-2 matrix in the form [C11 C12 ; 0 C22] where
% C11 is upper triangular (or "squeezed" upper trapezoidal), of dimension
% m1-by-n1.  The columns of C11 are the column singletons of A.  If C11 is
% square then all of its diagonal entries are larger in magnitude than tol.
%
% The input tol is optional; it defaults to 20*(m+n)*eps*maxcol2norm where
% [m n] = size(A) and maxcol2norm is the maximum 2-norm of the columns of A.
% The output tol is the tolerance used.
%
% If C11 is rectangular, then some of the column singletons have no
% corresponding row singleton.  A column j in C11 has a corresponding row i if
% if i = max(find(C(:,j))) > max(find(C(:,j-1))).  If present, abs(C(i,j)) will
% be larger than tol.
%
% Example:
%
%   A = [1 8 0 0
%        2 9 3 7
%        0 4 0 0
%        0 6 0 0 ]
%   [p q m1 n1 tol] = spqr_singletons (sparse (A))
%   C = A(p,q)
%
% In this example, C11 is 2-by-3.  One of the 3 column singletons (the 2nd on
% in C) has no corresponding row singleton.  This is an auxiliary routine that
% illustrates the singleton removal step used by SuiteSparseQR.  It is not need
% to solve a least-squares problem or find a sparse QR factorization; use
% SPQR_SOLVE or SPQR for those tasks, respectively.
%
% See also DMPERM, SPQR, SPQR_SOLVE.

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

error ('spqr_singletons mexFunction not found') ;


text:
classdef(Abstract) AbstractDragCoefficientModel < matlab.mixin.SetGet & matlab.mixin.Heterogeneous & matlab.mixin.Copyable
    %AbstractDragCoefficientModel Summary of this class goes here
    %   Detailed explanation goes here

    properties(Abstract, Constant)
        enum DragCoefficientModelEnum
    end

    methods
        CdA = getDragCoeff(obj, ut, rVect, vVect, bodyInfo, mass, altitude, pressure, density, vVectECEFMag, totalAoA, aoa, sideslip);

        tf = usesTotalAoA(obj)

        tf = usesAoaAndSideslip(obj)

        useTf = openEditDialog(obj)
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function installRes = plugin_askinstall(pluginName, pluginFunc, forceInstall)

if nargin < 3, forceInstall = false; end;

if nargin < 2 || ~exist(pluginFunc)

    if ~forceInstall
        db = dbstack;
        if length(db) > 2 && ~strcmpi(db(end).name, 'checkouteeglab.m');
            error([ 'Cannot find ' pluginName ' extension, use EEGLAB Extension Manager to install it' ]);
        end;

        installRes = 0;

        % check is deactivated
        try, PLUGINLIST = evalin('base', 'PLUGINLIST'); catch, PLUGINLIST = []; end;
        if ~isempty(PLUGINLIST) && isfield(PLUGINLIST, 'plugin')
            indPlugin = strmatch(lower(pluginName), lower({ PLUGINLIST.plugin }), 'exact');
            if ~isempty(indPlugin) && strcmpi(PLUGINLIST(indPlugin(1)).status, 'deactivated')
                res = questdlg2( [ pluginName ' extension is de-activated. Do you want to reactivate it now?' ], [ pluginName ' extension installation' ], 'No', 'Yes', 'Yes' );
                if strcmpi(res, 'no'), return, end;
                plugin_reactivate(PLUGINLIST(indPlugin(1)).foldername);
                evalin('base', 'eeglab rebuild');
                installRes = 1;
                return;
            end;
        end;

        % check for installing
        res = questdlg2( [ pluginName ' extension is not installed. Do you want to download it now?' ], [ pluginName ' extension installation' ], 'No', 'Yes', 'Yes' );
    else
        res = 'yes';
    end;

    if strcmpi(res, 'no'), return, end;
    plugins = plugin_getweb('import', []);
    indPlugin = strmatch(lower(pluginName), lower({ plugins.name }), 'exact');
    if isempty(indPlugin),
        plugins = plugin_getweb('process', []);
        indPlugin = strmatch(lower(pluginName), lower({ plugins.name }));
        if isempty(indPlugin),
            error([ pluginName ' extension not found' ]);
        end;
    end;
    result = plugin_install(plugins(indPlugin(1)).zip, plugins(indPlugin(1)).name, plugins(indPlugin(1)).version, forceInstall);
    if result == 1, installRes = 1; end;

    evalin('base', 'eeglab rebuild');
else
    installRes = 1;
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function F = in_fread_compumedics_pfs(sFile, sfid, SamplesBounds, ChannelsRange)
% IN_FREAD_COMPUMEDICS_PFS:  Read a block of recordings from a Compumedics ProFusion Sleep 4 exported binary file (.sdy/.rda).
%
% USAGE:  F = in_fread_compumedics_pfs(sFile, sfid, SamplesBounds=[], ChannelsRange=[])

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2015

% Parse inputs
if (nargin < 4) || isempty(ChannelsRange)
    ChannelsRange = [1, sFile.header.nchannels];
end
if (nargin < 3) || isempty(SamplesBounds)
    SamplesBounds = round(sFile.prop.times .* sFile.prop.sfreq);
end

% ===== COMPUTE OFFSETS =====
nChannels     = double(sFile.header.nchannels);
nReadTimes    = SamplesBounds(2) - SamplesBounds(1) + 1;
nReadChannels = double(ChannelsRange(2) - ChannelsRange(1) + 1);
% Everything is stored on 32 bit floats
bytesPerVal = 4;
dataClass = 'single';
% Header offset
offsetHeader = sFile.header.rda.segment(1).pos;
% Time offset
offsetTime = round(SamplesBounds(1) * nChannels * bytesPerVal);
% Channel offset at the beginning and end of each channel block
offsetChannelStart = round((ChannelsRange(1)-1) * bytesPerVal);
offsetChannelEnd   = (nChannels - ChannelsRange(2)) * bytesPerVal;
% Start reading at this point
offsetStart = offsetHeader + offsetTime + offsetChannelStart;
% Number of time samples to skip after each channel
offsetSkip = offsetChannelStart + offsetChannelEnd;

% ===== READ DATA BLOCK =====
% Position file at the beginning of the trial
fseek(sfid, offsetStart, 'bof');
% Read trial data
% => WARNING: CALL TO FREAD WITH SKIP=0 DOES NOT WORK PROPERLY
if (offsetSkip == 0)
    F = fread(sfid, [nReadChannels, nReadTimes], dataClass);
else
    precision = sprintf('%d*%s', nReadChannels, dataClass);
    F = fread(sfid, [nReadChannels, nReadTimes], precision, offsetSkip);
end
% Check that data block was fully read
if (numel(F) < nReadTimes * nReadChannels)
    % Error message
    disp(sprintf('BST> ERROR: File is truncated (%d values were read instead of %d)...', numel(F), nReadTimes * nReadChannels));
    % Pad with zeros
    Ftmp = zeros(nReadChannels, nReadTimes);
    Ftmp(1:numel(F)) = F(:);
    F = Ftmp;
end

author: brainstorm-tools


text:
function [res, bytesReq] = hlp_haveSufficientMemory(numelts,datatype)
% Check if we have sufficient memory to allocate an array.
%
% Inputs:
%   numelts: The number of elements in the array.
%       If the array is sparse, then numelts is the number of non-zero elements
%   datatype: The class of the array.
%       This can be any of the following: logical,int*,uint*,single,double,char
%
% Output:
%   res: true if we have sufficient memory, else false
%   bytesReq: the number of bytesRequired to allocate the array
%
% Author: Tim Mullen, SCCN/INC/UCSD, Jan, 2014

switch datatype
    case {'logical','int8','uint8'} , N=1;
    case {'char','int16','uint16'}  , N=2;
    case {'double','int64','uint64'}, N=8;
    case {'single','int32','uint32'}, N=4;
    otherwise
        error('unsupported datatype %s',datatype);
end

bytesReq = numelts*N;
res = (bytesReq <= hlp_getAvailableMemory('Bytes'));

repo: REST


text:
% example_radar.m
% ====================================================>
% The following seem to be descent parameters:
%   - ProbOfDetection = 0.8
%   - ProbOfDeath = 0.2
%   - VelocityErrVariance = 2
%   - ObsErrVariance = 50
%   - ProbOfConfirm = 0.9
%   - PHD.BirthScheme = {'Expansion', 5000}
%

%% Radar specific settings
RadarName = 'staddon';
switch(RadarName)
    case('staddon')
        % Surveillance region parameters
        V_bounds = [-8154.72944624983;... % X-min |
                    -212.289393440959;... % X-max | Surveillance region bounding
                    -7548.44272179096;... % Y-min | box coordinates (m)
                    4355.32645897434]';   % Y-max |
%         RadarCoords.lat = 50.33933333;
%         RadarCoords.lon = -4.12527778;
        RadarCoords.lat = 50.346069;
        RadarCoords.lon = -4.113670;
    case('longroom')
        % Surveillance region parameters
        V_bounds = [-8154.72944624983;... % X-min |
                    3000;...              % X-max | Surveillance region bounding
                    -7548.44272179096;... % Y-min | box coordinates (m)
                    4355.32645897434]';   % Y-max |
        RadarCoords.lat = 50.36286670714617;
        RadarCoords.lon = -4.156833300366998;
end

% Load dataset
load(strcat(RadarName,'.mat'));
N = size(DataList,2); % Simulation length

%% Focus on region of high clutter
% lon_lim = [-4.195 -4.11];
% lat_lim = [50.31 50.375];
% lon_lim = [-4.15 -4.13];
% lat_lim = [50.335 50.35];
% lon_lim = [-4.1538   -4.1152];
% lat_lim = [50.3097   50.35];

% lon_lim = [-4.1538   -4.1152];
% lat_lim = [50.3097   50.3373];

% lon_lim = [-4.167372 -4.148972];
% lat_lim = [50.353786, 50.361775];

lon_lim = [-4.1748   -4.1151];
lat_lim = [50.2992   50.3351];

% Convert LLA limits to NED
[y_lim, x_lim] = geodetic2ned(lat_lim,...
                              lon_lim, ...
                              0,...
                              RadarCoords.lat,...
                              RadarCoords.lon,...
                              0,...
                              referenceEllipsoid('wgs84'));
V_bounds = [x_lim(1);... % X-min |
            x_lim(2);... % X-max | Surveillance region bounding
            y_lim(1);... % Y-min | box coordinates (m)
            y_lim(2)]';  % Y-max |
for i = 1:N
    j = 1;
    while j <= size(DataList{i},2)
        if(DataList{i}(3,j)>lat_lim(2) || DataList{i}(3,j)<lat_lim(1) || DataList{i}(4,j)<lon_lim(1) || DataList{i}(4,j)>lon_lim(2))
            DataList{i}(:,j) = [];
            j = j -1;
        end
        j = j + 1;
    end
end

radar_meas_convert;
%% Plot & Recording settings
% Plot settings
ShowPlots = 1;              % Set to 0 to prevent showing any plots
ShowUpdate = 1;             % Set to 0 to skip showing update plots
ShowTrackInfo = 0;
NumPersistFrames = 10000;
PlotLatLon = 1;

% Recording settings
Record = 1;                 % Set to (0|1) to turn video recording (off|on)
FrameRate = 0.5;            % Number of frames per second
VideoQuality = 100;         % Set to desired quality percentage
VideoPathName = strcat(RadarName,'_heavy_clutter_staddon.avi'); % Set to the desired path and name of produced recording

% Model parameter shortcuts
lambdaV = 2; % Expected number of clutter measurements over entire surveillance region
V = (abs(V_bounds(2)-V_bounds(1))*abs(V_bounds(4)-V_bounds(3))); % Total area of surveillance region
P_D = 0.5;    % Probability of detection
timestep_duration = duration(0,0,2);

%% Instantiation of necessary components

% Instantiate a Dynamic model (CV model with q = 1 m/s^2)
% transition_model = ConstantVelocityX('VelocityErrVariance', 0.1,...
%                                      'NumDims', 2,...
%                                      'TimestepDuration', timestep_duration);
transition_model = OrnsteinUhlenbeckModelX('NumDims',2,...
                             'VelocityErrVariance',0.1,...
                             'DampingCoefficient',0.01,...
                             'TimestepDuration',timestep_duration);
% Instantiate an Observation model (Variance of 50m^2 on each coordinate)
measurement_model = LinearGaussianX('NumMeasDims', 2,...
                                    'NumStateDims', 4,...
                                    'MeasurementErrVariance', 10^2,...
                                    'Mapping', [1 3]);
clutter_model = PoissonRateUniformPosition2X('ClutterRate',lambdaV,...
                                            'Limits',[V_bounds(1:2);...
                                                      V_bounds(3:4)]);
detection_model = ConstantDetectionProbabilityX('DetectionProbability',P_D);

birth_model = DistributionBasedBirthModelX('Distribution', UniformDistributionX([V_bounds(1:2); ...
                                                                                [-7 7 ];...
                                                                                V_bounds(3:4);...
                                                                                [-7 7 ]]),...
                                           'BirthIntensity', 1);

% Compile the State-Space model
model = StateSpaceModelX(transition_model, measurement_model,...
                        'Clutter',clutter_model,...
                        'Detection', detection_model,...
                        'Birth', birth_model);

%% Base Filter
obs_covar= measurement_model.covar();
PriorState = GaussianStateX(zeros(4,1), transition_model.covar() + blkdiag(obs_covar(1,1), 0, obs_covar(2,2), 0));
base_filter = ExtendedKalmanFilterX('Model', model, 'StatePrior', PriorState);

%% Data Associator
config.ClutterModel = clutter_model;
config.Clusterer = NaiveClustererX();
config.Gater = EllipsoidalGaterX(2,'GateLevel',10)';
config.DetectionModel = detection_model;
jpdaf = JointIntegratedProbabilisticDataAssocX(config);

%% Track Initiator

% Initiate Data Associator
config_phd.Model = model;
[priorParticles, priorWeights] = model.Birth.random(50000);
config_phd.StatePrior = ParticleStateX(priorParticles,10*priorWeights);
config_phd.BirthScheme = {'Expansion', 5000};
config_phd.SurvivalProbability = 0.995;

% Instantiate PHD filter
myphd = SMC_PHDFilterX(config_phd);
%myphd = GM_PHDFilterX(config_phd2);
% Initiate Tag Generator
tag_gen = UuidTagGeneratorX();

% Prepare initiator parameters
config_ti.TagGenerator = tag_gen;
config_ti.InitFilter = base_filter;
config_ti.PhdFilter = myphd;
config_ti.ConfirmThreshold = 0.8;

% Create the track initiator
myti = PhdTrackInitiatorX(config_ti);

%% Track Deleter
config_td = struct('Fieldname', 'ExistenceProbability',...
                   'ReferenceValue', 0.1,...
                   'ReferenceOperand', 'lt');
mytd = FieldBasedDeleterX(config_td);

TrackList = [];

%% Create plot windows
if(ShowPlots)

    % Map plot
    figure('units','normalized','outerposition',[0 0 .5 0.9])
    ax(1) = gca;
    if PlotLatLon
        plot_google_map('Axis',ax(1),'APIKey','AIzaSyBXKujdtXRZiqya1soVS9pxBzYR4g7aGvM','Resize',3,'Scale',2,'MapType','satellite');
        axis(ax(1),[lon_lim(1) lon_lim(2) lat_lim(1) lat_lim(2)])
    end
%     axis(ax(1),[-4.195 -4.11 50.31 50.375])

%     % PHD Intensity plot
%     figure('units','normalized','outerposition',[.5 0 .5 1])
%     ax(2) = gca;

    plots = [];
end

%% START OF SIMULATION
% ===================>

for k=2:N
    fprintf('Iteration = %d/%d\n================>\n',k,N);

    %% Extract DataList at time k
    MeasurementList = MeasurementScans(k);
    timestamp_km1 = MeasurementScans(k-1).Timestamp;
    timestamp_k = MeasurementList.Timestamp;
    dt = timestamp_k - timestamp_km1;
    transition_model.TimestepDuration = dt;
    fprintf('Timestamp = %s\n================>\n',timestamp_k);

    %% Process JPDAF
    jpdaf.MeasurementList = MeasurementList;
    jpdaf.TrackList = TrackList;
    jpdaf.predictTracks();
    jpdaf.associate(TrackList, MeasurementList);
    jpdaf.updateTracks();

    % Delete all plots (other than map)
    for i = 1:numel(plots)
        delete(plots(i))
    end
    plots = [];
    hold on;

%     if(MeasurementList.NumMeasurements>0)
%         % Convert measurements to LLA and plot them
%         [lat,lon,~] = ned2geodetic(MeasurementList.Vectors(2,:),...
%                                    MeasurementList.Vectors(1,:),...
%                                    0,...
%                                    RadarCoords.lat,...
%                                    RadarCoords.lon,...
%                                    0,...
%                                    referenceEllipsoid('wgs84'));
%         plots(end+1) = plot(ax(1), lon,lat,'y*','MarkerSize', 10);
%
%         for i=1:numel(lat)
%             plots(end+1) = text(ax(1), lon(:,i)+0.0001,lat(:,i)+0.00027,num2str(i),'FontSize',8,'Color','w');
%         end
% %             plot(ax(1), RadarCoords.lon,RadarCoords.lat,...
% %                                '-s','MarkerSize',20,...
% %                                'MarkerEdgeColor','red',...
% %                                'MarkerFaceColor',[1 .6 .6]);
%     end
%     pause(0.01);
     %% Perform Track initiation
    [TrackList] = myti.initiateTracks(jpdaf.TrackList, MeasurementList, jpdaf.AssocWeightsMatrix);

    %% Perform Track deletion
    TrackList = mytd.deleteTracks(TrackList);

    % Plot update step results
    if(ShowPlots && ShowUpdate)

        % Delete all plots (other than map)
        for i = 1:numel(plots)
            delete(plots(i))
        end
        plots = [];
        hold on;

        if(MeasurementList.NumMeasurements>0)
            if PlotLatLon
                % Convert measurements to LLA and plot them
                [lat,lon,~] = ned2geodetic(MeasurementList.Vectors(2,:),...
                                           MeasurementList.Vectors(1,:),...
                                           0,...
                                           RadarCoords.lat,...
                                           RadarCoords.lon,...
                                           0,...
                                           referenceEllipsoid('wgs84'));
                plots(end+1) = plot(ax(1), lon,lat,'y*','MarkerSize', 10);
            else
                plots(end+1) = plot(ax(1), MeasurementList.Vectors(1,:),MeasurementList.Vectors(2,:),'b*','MarkerSize', 10);
            end
%             plot(ax(1), RadarCoords.lon,RadarCoords.lat,...
%                                '-s','MarkerSize',20,...
%                                'MarkerEdgeColor','red',...
%                                'MarkerFaceColor',[1 .6 .6]);

        end
        if PlotLatLon
            x = [-3560.53119790131, -3560.53119790131,  -3080.76025084906,  -3080.76025084906,  -3367.20269177091,  -3367.20269177091,  -3560.53119790131];
            y = [1349.54778996881,  1590.19881386064,   1590.19881386064,   1108.03444256592,   1108.03444256592,   1349.54778996881,   1349.54778996881];
            [y,x,~] = ned2geodetic(y,...
                                   x,...
                                   0,...
                                   RadarCoords.lat,...
                                   RadarCoords.lon,...
                                   0,...
                                   referenceEllipsoid('wgs84'));
            plots(end+1) = plot(ax(1), x, y, 'm-');
%             x = [-3367.20269177091,-3367.00101398088,-3083.23234936543,-3083.41702988195,-3367.20269177091];
%             y = [1108.20727096473,1425.35687149648,1425.18405344934,1108.03444256592,1108.20727096473];
%             [y,x,~] = ned2geodetic(y,...
%                                    x,...
%                                    0,...
%                                    RadarCoords.lat,...
%                                    RadarCoords.lon,...
%                                    0,...
%                                    referenceEllipsoid('wgs84'));
%             plots(end+1) = plot(ax(1), x, y, 'm-');
            x = [-2858.62556879160,-2856.97056906902,-108.925126527319,-108.988225201227,-2858.62556879160];
            y = [-4044.74840882411,-974.655039469575,-975.424226884065,-4045.51804181679,-4044.74840882411];
            [y,x,~] = ned2geodetic(y,...
                                   x,...
                                   0,...
                                   RadarCoords.lat,...
                                   RadarCoords.lon,...
                                   0,...
                                   referenceEllipsoid('wgs84'));
            plots(end+1) = plot(ax(1), x, y, 'r-');
        else
%             x = [-3560.53119790131, -3560.53119790131,  -3080.76025084906,  -3080.76025084906,  -3367.20269177091,  -3367.20269177091,  -3560.53119790131];
%             y = [1349.54778996881,  1590.19881386064,   1590.19881386064,   1108.03444256592,   1108.03444256592,   1349.54778996881,   1349.54778996881];
%             plots(end+1) = plot(ax(1), x, y, 'm-');
%             x = [-3367.20269177091,-3367.00101398088,-3083.23234936543,-3083.41702988195,-3367.20269177091];
%             y = [1108.20727096473,1425.35687149648,1425.18405344934,1108.03444256592,1108.20727096473];
%             plots(end+1) = plot(ax(1), x, y, 'm-');
            x = [-2858.62556879160,-2856.97056906902,-108.925126527319,-108.988225201227,-2858.62556879160];
            y = [-4044.74840882411,-974.655039469575,-975.424226884065,-4045.51804181679,-4044.74840882411];
            plots(end+1) = plot(ax(1), x, y, 'r-');
        end
        % Plot all existing tracks
        for j=1:numel(TrackList)
            track = TrackList{j};
            % Convert track trajectory to LLA and plot it
            means = [track.Trajectory.Mean];
            if PlotLatLon
                [lat,lon,~] = ned2geodetic(means(3,:),...
                                           means(1,:),...
                                           0,...
                                           RadarCoords.lat,...
                                           RadarCoords.lon,...
                                           0,...
                                           referenceEllipsoid('wgs84'));
                traj_length = size(lon,2);
                if(traj_length>NumPersistFrames)
                    start = traj_length-NumPersistFrames;
                else
                    start = 1;
                end
                plots(end+1) = plot(ax(1), lon(:,start:end),lat(:,start:end),'-.w','LineWidth',2);
                plots(end+1) = plot(ax(1), lon(:,end),lat(:,end),'ws','MarkerSize',15);
%                 plots(end+1) = gau

                % Convert track velocity to LLA and plot it
                [lat_vel,lon_vel,~] = ned2geodetic(track.State.Mean(4,end),...
                                                   track.State.Mean(2,end),...
                                                   0,...
                                                   lat(:,end),...
                                                   lon(:,end),...
                                                   0,...
                                                   referenceEllipsoid('wgs84'));
                lat_vel = lat_vel-lat(:,end);
                lon_vel = lon_vel-lon(:,end);
                plots(end+1) = quiver(ax(1), lon(:,end),lat(:,end),20*lon_vel,20*lat_vel,'r','LineWidth',1.5);
                if ShowTrackInfo
                    speed_kmph = sqrt(track.State.Mean(4,end)^2+track.State.Mean(2,end)^2)*3.6;
                    speed_knot = speed_kmph/1.852;
                    plots(end+1) = text(ax(1), lon(:,end)+0.001,lat(:,end)+0.00027,strcat("Sog:",num2str(speed_knot,2)," kt"),'FontSize',8,'Color','w');
                    plots(end+1) = text(ax(1), lon(:,end)+0.001,lat(:,end)-0.00027,strcat("PoE:",num2str(track.ExistenceProbability*100,3)," %"),'FontSize',8,'Color','w');
                end
            else
                traj_length = size(lon,2);
                if(traj_length>NumPersistFrames)
                    start = traj_length-NumPersistFrames;
                else
                    start = 1;
                end
                plots(end+1) = plot(ax(1), means(1,start:end),means(3,start:end),'-.k','LineWidth',2);
                plots(end+1) = plot(ax(1), means(1,end),means(3,end),'ks','MarkerSize',15);
%                 lat_vel = lat_vel-lat(:,end);
%                 lon_vel = lon_vel-lon(:,end);
%
                x_vel = track.State.Mean(3,end)*cos(track.State.Mean(4,end));
                y_vel = track.State.Mean(3,end)*sin(track.State.Mean(4,end));
                plots(end+1) = quiver(ax(1), means(1,end),means(2,end),20*x_vel,20*y_vel,'r','LineWidth',1.5);

                if ShowTrackInfo
                    speed_kmph = track.State.Mean(3,end)*3.6;
                    speed_knot = speed_kmph/1.852;
                    plots(end+1) = text(ax(1), means(1,end)+60,means(2,end)+50,strcat("Sog:",num2str(speed_knot,2)," kt"),'FontSize',8,'Color','k');
                    plots(end+1) = text(ax(1), means(1,end)+60,means(2,end)-50,strcat("PoE:",num2str(track.ExistenceProbability*100,3)," %"),'FontSize',8,'Color','k');
                end
            end

            % TODO: Convert heading and state covariance to LLA and plot them
            % [lat,lon,h] = ned2geodetic(North,East,0,50.346069,-4.113670,0,referenceEllipsoid('wgs84'));
            % h2 = plot_gaussian_ellipsoid(TrackList{j}.Filter.StateMean([1 3]), TrackList{j}.Filter.StateCovar([1 3],[1 3]),1,20,ax(1));
            % plots(end+1) = h2;
            % plots(end+1) = text(ax(1),TrackList{j}.Filter.StateMean(1)+20,TrackList{j}.Filter.StateMean(3)-5,int2str(TrackList{j}.TrackID));
            % plots(end+1) = text(ax(1),TrackList{j}.Filter.StateMean(1)+20,TrackList{j}.Filter.StateMean(3)-130,num2str(TrackList{j}.ProbOfExist,2));
        end

        % Add axis labels
%         xlabel('Longitude')
%         ylabel('Latitude')
set(gca,'visible','off')
        pause(0.0001)
        % Store video frame
        if(Record)
            F(k) = getframe(ax(1));
        end
    end

end

% Create video file and write to it
if(Record)
    F = F(2:end);
    vidObj = VideoWriter(char(VideoPathName));
    vidObj.Quality = VideoQuality;
    vidObj.FrameRate = FrameRate;
    open(vidObj);
    writeVideo(vidObj, F);
    close(vidObj);
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function value = r8vec_insignificant ( n, r, s )

%*****************************************************************************80
%
%% R8VEC_INSIGNIFICANT determines if an R8VEC is insignificant.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the vectors.
%
%    Input, real R(N), the vector to be compared against.
%
%    Input, real S(N), the vector to be compared.
%
%    Output, logical R8VEC_INSIGNIFICANT, is TRUE if S is insignificant
%    compared to R.
%
  value = 1;

  for i = 1 : n

    t = r(i) + s(i);
    tol = eps * abs ( r(i) );

    if ( tol < abs ( r(i) - t ) )
      value = 0;
      exit;
    end

  end

  return
end


text:
function dfdelgraphexclude
%DFDELGRAPHEXCLUDE Called when the exclusion dataset changes

% Copyright 2003-2004 The MathWorks, Inc.
% $Revision: 1.1.6.2 $  $Date: 2004/01/24 09:35:29 $

% Find the exclusion graph's figure window, and delete it
t = get(0,'ShowHiddenHandles');
set(0,'ShowHiddenHandles','on');
c = get(0,'Child');
f = findobj(c,'flat','Type','figure','Tag','dfexcludegraph');
set(0,'ShowHiddenHandles',t);
delete(f);

author: zouchuhang

repo: LayoutNet


text:
% calc_mse.m
% calculate mse error between two images

function E = calc_psnr(img1,img2);

[h w c] = size(img1);

diff = double(img1) - double(img2);
dobro = diff(:,:,1).*diff(:,:,1);
sum1 = sum(sum(dobro(:,:,1)));
e = sum1/(h*w);
E = 20*log10(255/sqrt(e));

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


function [sindex] = OrderingIndexes(varordering,varnames,newstrng)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 'OrderingIndexes' find the specified order of variables

% input:
% - varnames    = the list (names) of all the variables in the database
% - varordering = the names and the order of the variables in the VAR

% output:
% - sindex = the index in the varnames that correspond to the var ordering
% - index_XX = the index that the variable XXX has in the VAR

% Filippo Ferroni,
% Revised, 3/21/2018
% Revised, 9/11/2019
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<3
    newstrng='';
end

[~,sindex] = ismember(varordering,varnames);
for ii = 1 : size(varordering,2)
    check=0;
    for jjj = 1: size(varnames,2)
         if strcmp(deblank(varnames{jjj}),deblank(varordering{ii})) == 1,
             check = 1;
         end
    end
    if check == 1
        assignin('base', [newstrng 'index_' varordering{ii} ], ii);
    else
        warning(['I did not find ' varordering{ii} ' in varnames'])
    end
end


text:
classdef TestShapeTransformer
    %TestShapeTransformer

    methods (Static)
        function test_1
            img1 = cv.imread(fullfile(mexopencv.root(),'test','basketball1.png'), ...
                'Grayscale',true, 'ReduceScale',2);
            img2 = cv.imread(fullfile(mexopencv.root(),'test','basketball2.png'), ...
                'Grayscale',true, 'ReduceScale',2);

            detector = cv.AKAZE();
            [keypoints1, descriptors1] = detector.detectAndCompute(img1);
            [keypoints2, descriptors2] = detector.detectAndCompute(img2);
            matcher = cv.DescriptorMatcher('BFMatcher', ...
                'NormType',detector.defaultNorm());
            matches = matcher.match(descriptors1, descriptors2);
            pts1 = cat(1, keypoints1.pt);
            pts2 = cat(1, keypoints2.pt);

            for i=1:2
                if i==1
                    tps = cv.ShapeTransformer('AffineTransformer', ...
                        'FullAffine',true);
                    assert(isequal(tps.FullAffine, true));
                else
                    tps = cv.ShapeTransformer('ThinPlateSplineShapeTransformer', ...
                        'RegularizationParameter',25000);
                    assert(isequal(tps.RegularizationParameter, 25000));
                end

                tps.estimateTransformation(pts1, pts2, matches);

                [cost, pts3] = tps.applyTransformation(pts2);
                validateattributes(cost, {'numeric'}, {'scalar'});
                validateattributes(pts3, {'numeric'}, {'size',size(pts2)});

                img3 = tps.warpImage(img2);
                validateattributes(img3, {class(img2)}, {'size',size(img2)});
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [dx, dxGrad] = cstDyn(x,u)
% [dx, dxGrad] = cstDyn(x,u)
%
% Computes the dynamics (and gradients) for a 1d point mass on a
% friction-less plane with a force actuator.
%

% q = x(1,:);   %Position
dq = x(2,:);   %Velocity
ddq = u(1,:);  %Acceleration

dx = [dq;ddq];   %Pack up derivative of state

if nargout == 2   % Analytic gradients
    nTime = length(u);

    dqGrad = zeros(1,6,nTime); %6 = [time + pos + vel + force + slack];
    dqGrad(1,3,:) = 1; %gradient dq wrt dq

    ddqGrad = zeros(1,6,nTime);  %6 = [time + angle + rate +  force + slack];
    ddqGrad(1,4,:) = 1;  %gradient ddq wrt u

    dxGrad = cat(1, dqGrad, ddqGrad);

end

end

author: MatthewPeterKelly

repo: OptimTraj