text:
function [x,endPop,bPop,traceInfo] = ga(bounds,evalFN,evalOps,startPop,opts,...
termFN,termOps,selectFN,selectOps,xOverFNs,xOverOps,mutFNs,mutOps)
% GA run a genetic algorithm
% function [x,endPop,bPop,traceInfo]=ga(bounds,evalFN,evalOps,startPop,opts,
%                                       termFN,termOps,selectFN,selectOps,
%                                       xOverFNs,xOverOps,mutFNs,mutOps)
%
% Output Arguments:
%   x            - the best solution found during the course of the run
%   endPop       - the final population
%   bPop         - a trace of the best population
%   traceInfo    - a matrix of best and means of the ga for each generation
%
% Input Arguments:
%   bounds       - a matrix of upper and lower bounds on the variables
%   evalFN       - the name of the evaluation .m function
%   evalOps      - options to pass to the evaluation function ([NULL])
%   startPop     - a matrix of solutions that can be initialized
%                  from initialize.m
%   opts         - [epsilon prob_ops display] change required to consider two
%                  solutions different, prob_ops 0 if you want to apply the
%                  genetic operators probabilisticly to each solution, 1 if
%                  you are supplying a deterministic number of operator
%                  applications and display is 1 to output progress 0 for
%                  quiet. ([1e-6 1 0])
%   termFN       - name of the .m termination function (['maxGenTerm'])
%   termOps      - options string to be passed to the termination function
%                  ([100]).
%   selectFN     - name of the .m selection function (['normGeomSelect'])
%   selectOpts   - options string to be passed to select after
%                  select(pop,#,opts) ([0.08])
%   xOverFNS     - a string containing blank seperated names of Xover.m
%                  files (['arithXover heuristicXover simpleXover'])
%   xOverOps     - A matrix of options to pass to Xover.m files with the
%                  first column being the number of that xOver to perform
%                  similiarly for mutation ([2 0;2 3;2 0])
%   mutFNs       - a string containing blank seperated names of mutation.m
%                  files (['boundaryMutation multiNonUnifMutation ...
%                           nonUnifMutation unifMutation'])
%   mutOps       - A matrix of options to pass to Xover.m files with the
%                  first column being the number of that xOver to perform
%                  similiarly for mutation ([4 0 0;6 100 3;4 100 3;4 0 0])

% Binary and Real-Valued Simulation Evolution for Matlab
% Copyright (C) 1996 C.R. Houck, J.A. Joines, M.G. Kay
%
% C.R. Houck, J.Joines, and M.Kay. A genetic algorithm for function
% optimization: A Matlab implementation. ACM Transactions on Mathmatical
% Software, Submitted 1996.
%
% 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 1, or (at your option)
% any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. A copy of the GNU
% General Public License can be obtained from the
% Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

%%$Log: ga.m,v $
%Revision 1.10  1996/02/02  15:03:00  jjoine
% Fixed the ordering of imput arguments in the comments to match
% the actual order in the ga function.
%
%Revision 1.9  1995/08/28  20:01:07  chouck
% Updated initialization parameters, updated mutation parameters to reflect
% b being the third option to the nonuniform mutations
%
%Revision 1.8  1995/08/10  12:59:49  jjoine
%Started Logfile to keep track of revisions
%

n=nargin;
if n<2 | n==6 | n==10 | n==12
  disp('Insufficient arguements')
end
if n<3 %Default evalation opts.
  evalOps=[];
end
if n<5
  opts = [1e-6 1 0];
end
if isempty(opts)
  opts = [1e-6 1 0];
end

if any(evalFN<48) %Not using a .m file
  if opts(2)==1 %Float ga
    e1str=['x=c1; c1(xZomeLength)=', evalFN ';'];
    e2str=['x=c2; c2(xZomeLength)=', evalFN ';'];
  else %Binary ga
    e1str=['x=b2f(endPop(j,:),bounds,bits); endPop(j,xZomeLength)=',...
	evalFN ';'];
  end
else %Are using a .m file
  if opts(2)==1 %Float ga
    e1str=['[c1 c1(xZomeLength)]=' evalFN '(c1,[gen evalOps]);'];
    e2str=['[c2 c2(xZomeLength)]=' evalFN '(c2,[gen evalOps]);'];
  else %Binary ga
    e1str=['x=b2f(endPop(j,:),bounds,bits);[x v]=' evalFN ...
	'(x,[gen evalOps]); endPop(j,:)=[f2b(x,bounds,bits) v];'];
  end
end

if n<6 %Default termination information
  termOps=[100];
  termFN='maxGenTerm';
end
if n<12 %Default muatation information
  if opts(2)==1 %Float GA
  mutFNs=['boundaryMutation multiNonUnifMutation nonUnifMutation unifMutation'];
    mutOps=[4 0 0;6 termOps(1) 3;4 termOps(1) 3;4 0 0];
  else %Binary GA
    mutFNs=['binaryMutation'];
    mutOps=[0.05];
  end
end
if n<10 %Default crossover information
  if opts(2)==1 %Float GA
    xOverFNs=['arithXover heuristicXover simpleXover'];
    xOverOps=[2 0;2 3;2 0];
  else %Binary GA
    xOverFNs=['simpleXover'];
    xOverOps=[0.6];
  end
end
if n<9 %Default select opts only i.e. roullete wheel.
  selectOps=[];
end
if n<8 %Default select info
  selectFN=['normGeomSelect'];
  selectOps=[0.08];
end
if n<6 %Default termination information
  termOps=[100];
  termFN='maxGenTerm';
end
if n<4 %No starting population passed given
  startPop=[];
end
if isempty(startPop) %Generate a population at random
  %startPop=zeros(80,size(bounds,1)+1);
  startPop=initializega(80,bounds,evalFN,evalOps,opts(1:2));
end

if opts(2)==0 %binary
  bits=calcbits(bounds,opts(1));
end

xOverFNs=parse(xOverFNs);
mutFNs=parse(mutFNs);

xZomeLength  = size(startPop,2); 	%Length of the xzome=numVars+fittness
numVar       = xZomeLength-1; 		%Number of variables
popSize      = size(startPop,1); 	%Number of individuals in the pop
endPop       = zeros(popSize,xZomeLength); %A secondary population matrix
c1           = zeros(1,xZomeLength); 	%An individual
c2           = zeros(1,xZomeLength); 	%An individual
numXOvers    = size(xOverFNs,1); 	%Number of Crossover operators
numMuts      = size(mutFNs,1); 		%Number of Mutation operators
epsilon      = opts(1);                 %Threshold for two fittness to differ
oval         = max(startPop(:,xZomeLength)); %Best value in start pop
bFoundIn     = 1; 			%Number of times best has changed
done         = 0;                       %Done with simulated evolution
gen          = 1; 			%Current Generation Number
collectTrace = (nargout>3); 		%Should we collect info every gen
floatGA      = opts(2)==1;              %Probabilistic application of ops
display      = opts(3);                 %Display progress

while(~done)
  %Elitist Model
  [bval,bindx] = max(startPop(:,xZomeLength)); %Best of current pop
  best =  startPop(bindx,:);

  if collectTrace
    traceInfo(gen,1)=gen; 		          %current generation
    traceInfo(gen,2)=startPop(bindx,xZomeLength);       %Best fittness
    traceInfo(gen,3)=mean(startPop(:,xZomeLength));     %Avg fittness
    traceInfo(gen,4)=std(startPop(:,xZomeLength));
  end

  if ( (abs(bval - oval)>epsilon) | (gen==1)) %If we have a new best sol
    if display
      fprintf(1,'\n%d %f\n',gen,bval);          %Update the display
    end
    if floatGA
      bPop(bFoundIn,:)=[gen startPop(bindx,:)]; %Update bPop Matrix
    else
      bPop(bFoundIn,:)=[gen b2f(startPop(bindx,1:numVar),bounds,bits)...
	  startPop(bindx,xZomeLength)];
    end
    bFoundIn=bFoundIn+1;                      %Update number of changes
    oval=bval;                                %Update the best val
  else
    if display
      fprintf(1,'%d ',gen);	              %Otherwise just update num gen
    end
  end

  endPop = feval(selectFN,startPop,[gen selectOps]); %Select

  if floatGA %Running with the model where the parameters are numbers of ops
    for i=1:numXOvers,
      for j=1:xOverOps(i,1),
	a = round(rand*(popSize-1)+1); 	%Pick a parent
	b = round(rand*(popSize-1)+1); 	%Pick another parent
	xN=deblank(xOverFNs(i,:)); 	%Get the name of crossover function
	[c1 c2] = feval(xN,endPop(a,:),endPop(b,:),bounds,[gen xOverOps(i,:)]);

	if c1(1:numVar)==endPop(a,(1:numVar)) %Make sure we created a new
	  c1(xZomeLength)=endPop(a,xZomeLength); %solution before evaluating
	elseif c1(1:numVar)==endPop(b,(1:numVar))
	  c1(xZomeLength)=endPop(b,xZomeLength);
	else
	  %[c1(xZomeLength) c1] = feval(evalFN,c1,[gen evalOps]);
	  eval(e1str);
	end
	if c2(1:numVar)==endPop(a,(1:numVar))
	  c2(xZomeLength)=endPop(a,xZomeLength);
	elseif c2(1:numVar)==endPop(b,(1:numVar))
	  c2(xZomeLength)=endPop(b,xZomeLength);
	else
	  %[c2(xZomeLength) c2] = feval(evalFN,c2,[gen evalOps]);
	  eval(e2str);
	end

	endPop(a,:)=c1;
	endPop(b,:)=c2;
      end
    end

    for i=1:numMuts,
      for j=1:mutOps(i,1),
	a = round(rand*(popSize-1)+1);
	c1 = feval(deblank(mutFNs(i,:)),endPop(a,:),bounds,[gen mutOps(i,:)]);
	if c1(1:numVar)==endPop(a,(1:numVar))
	  c1(xZomeLength)=endPop(a,xZomeLength);
	else
	  %[c1(xZomeLength) c1] = feval(evalFN,c1,[gen evalOps]);
	  eval(e1str);
	end
	endPop(a,:)=c1;
      end
    end

  else %We are running a probabilistic model of genetic operators
    for i=1:numXOvers,
      xN=deblank(xOverFNs(i,:)); 	%Get the name of crossover function
      cp=find(rand(popSize,1)<xOverOps(i,1)==1);
      if rem(size(cp,1),2) cp=cp(1:(size(cp,1)-1)); end
      cp=reshape(cp,size(cp,1)/2,2);
      for j=1:size(cp,1)
	a=cp(j,1); b=cp(j,2);
	[endPop(a,:) endPop(b,:)] = feval(xN,endPop(a,:),endPop(b,:),...
	  bounds,[gen xOverOps(i,:)]);
      end
    end
    for i=1:numMuts
      mN=deblank(mutFNs(i,:));
      for j=1:popSize
	endPop(j,:) = feval(mN,endPop(j,:),bounds,[gen mutOps(i,:)]);
	eval(e1str);
      end
    end
  end

  gen=gen+1;
  done=feval(termFN,[gen termOps],bPop,endPop); %See if the ga is done
  startPop=endPop; 			%Swap the populations

  [bval,bindx] = min(startPop(:,xZomeLength)); %Keep the best solution
  startPop(bindx,:) = best; 		%replace it with the worst
end

[bval,bindx] = max(startPop(:,xZomeLength));
if display
  fprintf(1,'\n%d %f\n',gen,bval);
end

x=startPop(bindx,:);
if opts(2)==0 %binary
  x=b2f(x,bounds,bits);
  bPop(bFoundIn,:)=[gen b2f(startPop(bindx,1:numVar),bounds,bits)...
      startPop(bindx,xZomeLength)];
else
  bPop(bFoundIn,:)=[gen startPop(bindx,:)];
end

if collectTrace
  traceInfo(gen,1)=gen; 		%current generation
  traceInfo(gen,2)=startPop(bindx,xZomeLength); %Best fittness
  traceInfo(gen,3)=mean(startPop(:,xZomeLength)); %Avg fittness
end

author: HuangCongQing


text:
function pwl_interp_2d_scattered_test01 ( )

%*****************************************************************************80
%
%% PWL_INTERP_2D_TEST01 tests R8TRIS2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 October 2012
%
%  Author:
%
%    John Burkardt
%
  node_num = 9;
  element_order = 3;

  node_xy = [ ...
       0.0, 0.0;
       0.0, 1.0;
       0.2, 0.5;
       0.3, 0.6;
       0.4, 0.5;
       0.6, 0.4;
       0.6, 0.5;
       1.0, 0.0;
       1.0, 1.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PWL_INTERP_2D_TEST01\n' );
  fprintf ( 1, '  R8TRIS2 computes the Delaunay triangulation of\n' );
  fprintf ( 1, '  a set of nodes in 2D.\n' );
%
%  Set up the Delaunay triangulation.
%
  [ element_num, triangle, element_neighbor ] = r8tris2 ( node_num, node_xy );

  triangulation_order3_print ( node_num, element_num, node_xy, ...
    triangle, element_neighbor );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function val=secDerivVal(z,n)
%%SECDERIVVAL Return the value of the nth derivative of sec(x) with respect
%             to x evaluated at x=z.
%
%INPUTS: z A matrix of real or complex values at which the nth derivative
%          of the secant function is desired.
%        n The number of derivatives to take. n>=0.
%
%OUTPUTS: val The value of the nth derivative of the secant function taken
%             at all of the points in z. val has the same dimensions as z.
%
%This function implements the algorithm of [1].
%
%REFERENCES:
%[1] M. E. Hoffman, "Derivative polynomials for tangent and secant," The
%    American Mathematical Monthly, vol. 102, no. 1, pp. 23-30, Jan. 1995.
%
%October 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numEls=numel(z);

val=zeros(size(z));
for curEl=1:numEls
    x=z(curEl);
    u=tan(x);

    P=zeros(n+1,1);
    Q=zeros(n+1,1);
    P(1)=u;
    Q(1)=1;

    for k=1:n
        %The loop implements Equation 5.
        sumValP=0;
        sumValQ=0;
        for i=0:(k-1)
            binomVal=binomial(k-1,i);
            sumValP=sumValP+binomVal*P(i+1)*P(k-1-i+1);
            sumValQ=sumValQ+binomVal*P(i+1)*Q(k-1-i+1);
        end
        if(k==1)
            sumValP=sumValP+1;
        end
        P(k+1)=sumValP;
        Q(k+1)=sumValQ;
    end
    val(curEl)=Q(end)*sec(x);
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [ n_data, mu, sigma, a, b, x, fx ] = ...
  truncated_normal_ab_pdf_values ( n_data )

%*****************************************************************************80
%
%% TRUNCATED_NORMAL_AB_PDF_VALUES: values of the Truncated Normal AB PDF.
%
%  Discussion:
%
%    The Normal distribution, with mean Mu and standard deviation Sigma,
%    is truncated to the interval [A,B].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real MU, the mean of the distribution.
%
%    Output, real SIGMA, the standard deviation of the distribution.
%
%    Output, real A, B, the lower and upper truncation limits.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 11;

  a_vec = [ ...
     50.0, ...
     50.0, ...
     50.0, ...
     50.0, ...
     50.0, ...
     50.0, ...
     50.0, ...
     50.0, ...
     50.0, ...
     50.0, ...
     50.0 ];

  b_vec = [ ...
     150.0, ...
     150.0, ...
     150.0, ...
     150.0, ...
     150.0, ...
     150.0, ...
     150.0, ...
     150.0, ...
     150.0, ...
     150.0, ...
     150.0 ];

  fx_vec = [ ...
    0.01543301171801836, ...
    0.01588394472270638, ...
    0.01624375997031919, ...
    0.01650575046469259, ...
    0.01666496869385951, ...
    0.01671838200940538, ...
    0.01666496869385951, ...
    0.01650575046469259, ...
    0.01624375997031919, ...
    0.01588394472270638, ...
    0.01543301171801836 ];

  mu_vec = [ ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0, ...
     100.0 ];

  sigma_vec = [ ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0, ...
    25.0  ];

  x_vec = [ ...
     90.0, ...
     92.0, ...
     94.0, ...
     96.0, ...
     98.0, ...
    100.0, ...
    102.0, ...
    104.0, ...
    106.0, ...
    108.0, ...
    110.0 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    a = 0.0;
    b = 0.0;
    mu = 0.0;
    sigma = 0.0;
    x = 0.0;
    fx = 0.0;
  else
    a = a_vec(n_data);
    b = b_vec(n_data);
    mu = mu_vec(n_data);
    sigma = sigma_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function [] = sb_fix_unwrap_manual(reset_flag,ix_reduced_list)
% [] = sb_fix_unwrap_manual(reset_flag)
% script to manually fix unwrapping errors in Small Baseline processing
%
% This function allows to manually fix unwrapping errors by iterating
% the results of ps_plot('rsb'). You can add add or subtract an integer number
% of 2pi cycles over a selected region, or correct with respect to the closest
% wrap of another region. You can keep iterating till the 'rsb'
% errors have reduced. If it goes wrong and you do want to reset your
% unwrapped data back to that of the original you can set the reset_flag to 1;
%
% By Bekaert David - University of Leeds - Sept 2014
%
% modifications:
% Bekaert David     05/2015     Fix in command line output, fix the code for ps_plot
% Bekaert David     09/2015     Add file restoring.
% Bekaert David     12/2016     Add closest wrap option, change to stamps_save code
% Bekaert David     12/2016     Add argument which is a subset of intergerograms to reduce plotting needs.
%

% deramp IFGs
deramp_flag =0;         % optional deramp rsb. Might be needed when processed with gamma
% what to use as guidence 'rsb' or 'usb' to click on
plot_option = 'usb';    % only 'rsb' or 'usb'
dem_option = 'd';    % only 'd' or ''

% flags
if nargin<1
    reset_flag=0;
end
if isempty(reset_flag)
    reset_flag=0;
end
if nargin<2
    ix_reduced_list=[];
end
if ~strcmpi(plot_option,'usb') && ~strcmpi(plot_option,'rsb')
    error('Only plot_option: rsb and usb supported');
end
if ~strcmpi(dem_option,'d') && ~isempty(dem_option)
    error('Only dem_option: d and [] supported');
end

% keeping the original data, incase user want to reset it
if exist('phuw_sb2_original.mat','file')~=2
    copyfile('phuw_sb2.mat','phuw_sb2_original.mat')
    copyfile('phuw2.mat','phuw2_original.mat')
    copyfile('phuw_sb_res2.mat','phuw_sb_res2_original.mat')
else
    copyfile('phuw_sb2.mat','phuw_sb2_previous_temp.mat')
    copyfile('phuw2.mat','phuw2_previous_temp.mat')
    copyfile('phuw_sb_res2.mat','phuw_sb_res2_previous_temp.mat')
end

% use previous or start from scratch
if reset_flag==1
   fprintf('Using the original data\n')
   ph_input = load('phuw_sb2_original');
else
   fprintf('Use data from a previous run\n')
   ph_input = load('phuw_sb2.mat');
end

% loading the interferogram information
ps = load('ps2.mat');

% generating list of ifgs to be ploted
drop_ifg_index = getparm('drop_ifg_index');
if ~isempty(ix_reduced_list)
    % removing ifgs that have been dropped before
    for k=1:length(drop_ifg_index)
        ix_temp = find(drop_ifg_index(k)==ix_reduced_list);
        ix_reduced_list(ix_temp)=[];
    end
    if sum(ix_reduced_list>ps.n_ifg)>0
        fprintf('Your list is larger than number of IFGS, will reset to max number of IFG \n')
        ix_reduced_list(ix_reduced_list>ps.n_ifg)=[];
    end
end
% reset the list in case nothing was left
if isempty(ix_reduced_list)
   ix_reduced_list = 1:ps.n_ifg;
   ix_reduced_list(drop_ifg_index)=[];
end

% deramping
if deramp_flag==1
    deramp_option = ['-o'];
else
    deramp_option = '';
end
if ~isempty(dem_option)
    if deramp_flag==1
        dem_option='-do';
    else
        dem_option= '-d';
    end
end
if strcmpi(plot_option,'usb')
    plot_option = [plot_option dem_option];
elseif strcmpi(plot_option,'rsb')
    plot_option = [plot_option deramp_option];
end

% plotting the current rsb data
ps_plot(['rsb' deramp_option],1,0,0,ix_reduced_list);

% get the interferogram that the user needs to adapt
repeat=1;
while repeat==1
    ix_ifg = input('Which interferogram to you want to correct? ','s');
    ix_ifg = str2num(ix_ifg);
    if isempty(ix_ifg)
        repeat=1;
    elseif ix_ifg<=ps.n_ifg
        repeat=0;
    else
        fprintf(['Not that many interferograms \n'])
    end
end

% plot the rsb value for this interferogram
% option one can use deramped rsb, this is for teh case teh interferograms
% were not created from relative differences, i.e. each interferogram had a
% baseline estimated and there might be some ramping errors because of that.
h_fig = ps_plot([plot_option],1,0,0,ix_ifg);
set(h_fig,'name',['Original interferogram']);

% Getting a polygon of the incorrect unwrapped area
fprintf('Define the incorrect unwrapped region through a polygon by clicking on the figure. \n')
repeat_zoom=1;
while repeat_zoom==1
    action_flag= input('First, zoom to your region of interest. Press [c] to continue. ','s');
    if strcmpi(action_flag,'c')
        repeat_zoom=0;
    end
end

fprintf('Now, start defining the polygon by outlining the incorrect region. \n Press enter once done\n')
% call the figure in case the user clicked somewhere else
figure(h_fig);
polygon=ginput;
% plotting the polygon on top
hold on
plot([polygon(:,1) ;polygon(1,1)],[polygon(:,2);polygon(1,2)],'r-','linewidth',2)

% loop untill the user is happy with it
continue_flag=1;
while continue_flag==1
    repeat=1;
    fprintf('You can shift the whole region by an integer number of cycles or you can put all pixels to a specifc wrap \n')
    while repeat==1
        ix_shift= input('By how many cycles to you want to shift this region? [+-integer or inf for wrap option] ','s');
        ix_shift = str2num(ix_shift);
        if isempty(ix_shift)
            repeat=1;
        elseif ix_shift==inf
            repeat=0;
        elseif (ix_shift./(round(ix_shift)))~=1
            fprintf(['Needs to be an integer number... \n'])
            repeat =1;
        else
            repeat=0;
        end
    end

    % finding the pixels within the polygon
    ix = inpolygon(ps.lonlat(:,1),ps.lonlat(:,2),polygon(:,1),polygon(:,2));

    % checking the option the user picked - closes wrap (inf) or shift region
    if ix_shift==inf            % closest wrap option
        fprintf('Define region to which you want to define as reference wrap (average will be used!). \n Press enter once done\n')
        % call the figure in case the user clicked somewhere else
        figure(h_fig);
        polygon_ref=ginput;
        % plotting the ploygon on top
        hold on
        plot([polygon_ref(:,1) ;polygon_ref(1,1)],[polygon_ref(:,2);polygon_ref(1,2)],'b-','linewidth',4)

        % finding the pixels within the reference polygon
        ix_ref = inpolygon(ps.lonlat(:,1),ps.lonlat(:,2),polygon_ref(:,1),polygon_ref(:,2));

        % check to which interferograms this should be applied
        repeat2=1;
        while repeat2==1
            action_flag= input('Do you want to apply this to all interferograms [y/n]? ','s');
            if strcmpi(action_flag,'y')
                repeat2=0;
            elseif strcmpi(action_flag,'n')
                repeat2=0;
            end
        end

        % store the orginal interferograms
        ix_ifg_or = ix_ifg;
        if strcmpi(action_flag,'y')
           ix_ifg=[1:size(ph_uw,2)];
        end

        % do the estimation for each itnerferogram
        ph_uw= ph_input.ph_uw;
        ref_phase = nanmean(ph_uw(ix_ref,ix_ifg),1);
        for k_ifgs=1:length(ix_ifg)
           % compute the reference
            radian_shift = round((ph_uw(ix,ix_ifg(k_ifgs))-ref_phase(k_ifgs))./(2*pi))*2*pi;
            ph_uw(ix,ix_ifg(k_ifgs)) = ph_uw(ix,ix_ifg(k_ifgs)) - radian_shift;
        end

        % update back to the previous interferogram to be corrected for
        % plotting purposes
        ix_ifg = ix_ifg_or;
        clear ref_phase

    else        % option of shifting the interferogram
        % the shift in radians
        radian_shift = ix_shift*2*pi;

        % modifying the interferogram
        ph_uw= ph_input.ph_uw;
        ph_uw(ix,ix_ifg)=ph_uw(ix,ix_ifg)+radian_shift;
    end
    msd = ph_input.msd;

    % saving the data
    stamps_save('phuw_sb2.mat',ph_uw,msd);

    % re-running the
    sb_invert_uw

    % plot the new residuals
    ps_plot(['rsb' deramp_option],1,0,0,ix_reduced_list);
    h_fig_new = ps_plot([plot_option],1,0,0,ix_ifg);
    set(h_fig_new,'name',['Corrected interferogram']);

    repeat=1;
    while repeat==1
        string= input('retry? [y/n] ','s');
        if strcmpi(string,'y')
            repeat=0;
            continue_flag = 1;
        elseif strcmpi(string,'n')
            repeat=0;
            continue_flag = 0;
            % see if the result needs to be kept or reverted
            repeat2=1;
            while repeat2==1
                action_flag= input('Keep this result [y/n]? ','s');
                if strcmpi(action_flag,'y')
                    repeat2=0;
                elseif strcmpi(action_flag,'n')
                    repeat2=0;

                    % restore the codes
                    copyfile('phuw_sb2_original.mat','phuw_sb2.mat')
                    copyfile('phuw2_original.mat','phuw2.mat')
                    copyfile('phuw_sb_res2_original.mat','phuw_sb_res2.mat')

                else
                    fprintf('y or n ...\n')
                end
            end
        else
            fprintf('y or n ...\n')
        end
    end

end

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function diffVal=angCircDiff(ang)
%%ANGCIRCDIFF Given two angles (in radians) on a unit circle, take their
%         difference in radians. The difference can only span +/-pi. For
%         example the difference between 0 and 2*pi should be zero as both
%         represent the same point.
%
%INPUTS: angs A 2XnumAng matrix of angles in radians. The differences
%               angs(1,:)-angs(2,:) in each column will be taken.
%
%OUTPUTS: diffVal A 1XnumAng vector of the differences on the unit circle
%                 of the angles. The shortest way around the circle is
%                 chosen. This goes from -pi to pi.
%
%This function simple takes the difference and wraps it to the range of -pi
%to pi with the wrapRange function.
%
%EXAMPLE:
%One can verify that for an angle of 0.1
% angCircDiff([0;0.1])-angCircDiff([0;0.1+2*pi])
%the difference is around machine precision. Similarly,
% angCircDiff([0;2*pi])
% is exactly zero
%
%January 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

diffVal=wrapRange(ang(1,:)-ang(2,:),-pi,pi);

end

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


text:
function state = viterbi(feature, hmm)
%VITERBI    Viterbi decoding of HMM.
%   Format: state = viterbi(feature, hmm)
%   Inputs:
%       feature:    Feature matrix, where each column stands for a time
%                   step.
%       HMM:        HMM model.
%   Output:
%       state:      Decoded state sequence.

    frames = size(feature, 2);
    states = length(hmm.init);
    decision = zeros(frames, states);

    emit = -inf(frames, states);
    for u = 1:states
        if (~isempty(hmm.gmm{u}))
            emit(:,u) = log(hmm.gmm{u}.pdf(feature'));
        end
    end

    trans = log(hmm.trans);
    prob = log(hmm.init);
    for f = 1:frames
        newProb = zeros(1, states);
        for u = 1:states
            [newProb(u) decision(f,u)] = max(prob + trans(:,u)');
        end
        prob = newProb + emit(f,:);
    end

    [temp state(frames)] = max(prob);
    for f = frames:-1:2
        state(f-1) = decision(f, state(f));
    end
end

author: MaigoAkisame


text:
function label_mat  = convert_labelvec_to_mat(label_vec, num, class_num)
    label_mat = zeros(class_num, num);
    for i=1:num
        label_mat(label_vec(i),i) = 1;
    end
end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
% Convert a stroke [x,y] such that it is uniformly sampled in space.
%
% Input
%   stk : [n x 2] stroke
%   dint : target distance between poitns
%
% Output
%   yi : [m x 2] interpolated stroke
%

%%
function stk_yi = uniform_space_lerp(stk,dint)

    %% return if stroke is too short
    n = size(stk,1);
    if n==1
       stk_yi = stk;
       return;
    end

    %% compute distance between each point
    dist = zeros(n,1);
    tormv = false(n,1);
    for i=2:n
        x1 = stk(i,:);
        x2 = stk(i-1,:);
        dist(i) = norm(x1-x2);
        tormv(i) = dist(i) < 1e-4;
    end

    %% remove points that are too close
    dist(tormv) = [];
    stk(tormv,:) = [];

    %% return if stroke is too short
    n = size(stk,1);
    if n==1
       stk_yi = stk;
       return;
    end

    %% cumulative distance
    cumdist = cumsum(dist);
    start_dist = cumdist(1);
    end_dist = cumdist(end);
    x = cumdist(:);

    %%
    nint = round(end_dist/dint);
    nint = max(nint,2);
    xi = linspace(start_dist,end_dist,nint);
    stk_yi = interp1(x,stk,xi);

    bool_viz = false;
    if bool_viz
       figure
       hold on
       plot(stk_yi(:,1),stk_yi(:,2),'b.','MarkerSize',18);
       plot(stk(:,1),stk(:,2),'r.','MarkerSize',6);
       xlim([0 105]);
       ylim([-105 0]);
    end

end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
%parademo.m
%
% REQUIRES THE DATA SETS !!!!!
%

% Copyright (C) 1995-2006  Rasmus Bro & Claus Andersson
% Copenhagen University, DK-1958 Frederiksberg, Denmark, rb@life.ku.dk
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 2 of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
% Street, Fifth Floor, Boston, MA  02110-1301, USA.

close all
clear all
echo off
home
echo on
load claus;
echo off

disp(' ')
disp(' Fluorescence measurements ideally follow the trilinear PARAFAC')
disp(' model. We will use a simple data set of 5 mixtures of three')
disp(' amino-acids (trp,phe & tyr) to show how the pure spectra and')
disp(' concentrations can be found with parafac')
disp(' ')
disp(' First lets look at the raw data:')
disp(' Press any key to continue')
pause

echo on
figure(1);
for i=1:5,
   subplot(3,2,i)
   sample = squeeze(X(i,:,:));
   mesh(ExAx,EmAx,sample);
   title(['Raw data - sample ',num2str(i)]);
   xlabel('Excitation [nm]')
   ylabel('Emission [nm]')
   axis([ExAx(1) ExAx(end) EmAx(1) EmAx(end) 0 1000]);
   grid on
   drawnow
end;
echo off

disp(' ')
disp(' Press any key to continue')
pause
close all
home

disp(' ')
disp(' The data will be slightly reduced to save computation time. This is done')
disp(' by reducing the 2. and 3. emission and excitation mode:')
disp(' ')

echo on
X = X(:,1:5:end,1:2:end);
EmAx = EmAx(1:5:end);
ExAx = ExAx(1:5:end);
size(X)
echo off

disp(' ')
disp(' Press any key to continue')
pause
close all
home

disp(' ')
disp(' PARAFAC may be fitted to the data, but in order to investigate how')
disp(' many components are needed we will use the tool pftest to get an')
disp(' indication. We use 1 to 5 components because we assume that appr.')
disp(' 3 are reasonable but want to check that. We fit models from 1 to 5 ')
disp(' components and fit each three times. This way we can check that the')
disp(' is the same every time we fit, say a four-component model. If it is')
disp(' it is an indication that we have used too many components or that the')
disp('  data are difficult to fit')
disp(' ')
disp(' This process will take some time, but afterwards a plot is produced')
disp(' which is often very instructive to look at')
disp(' ')
disp(' ')
disp(' Press any key to continue')
disp(' ')
pause
close all
home

echo on
[ssX,Corco] = pftest(3,X,5,[0 0 0 0 NaN]);
echo off
disp(' ')
disp(' Indeed, the fit values seem to indicate that three components are ')
disp(' suitable, whereas the core consistency seems to point to a possible')
disp(' fourth component. This fourth component must be weak though considering')
disp(' the low increase in fit. If looking into the 3- and 4-component models')
disp(' it is realized that the fourth component is reflecting Rayleigh scatter.')
disp(' The problem could have been avoided if we remove emission below ')
disp(' excitation (by setting these elements to NaN), but this is not pursued')
disp(' here. Instead, we will fit a three-component parafac model.')
disp(' ')
disp(' We will turn the plotting on, so that a graphical output is produced')
disp(' ')
disp(' Press any key to continue')
pause
close all
home

echo on
model = parafac(X,3,[0 0 2]);
echo off

disp(' ')
disp(' Press any key to continue')
pause
close all
home

disp(' Using the function fac2let, we can get scores and loading out')
disp(' ')
echo on
[A,B,C] = fac2let(model);
echo off
disp(' ')
disp(' Scores')
disp(A)
disp(' ')
disp(' Concentrations')
disp(y)
disp(' ')
disp(readme)
disp(' ')

disp(' Scrutinize the scores and compare to the concentrations to see')
disp(' the relation between the two. Each score correspond to a specific')
disp(' component up to scale and permutation')
disp(' ')
disp(' END OF  PARADEMO')

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
classdef PTKGetLeftLungROI < PTKPlugin
    % PTKGetLeftLungROI. Plugin for finding the left lung region of interest.
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %     PTKGetLeftLungROI runs the library function PTKLeftAndRightLungs to
    %     segment the left and right lungs, and uses this to determine a region
    %     of interest for the left lung only.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Left Lung <BR>ROI'
        ToolTip = 'Change the context to display the left lung'
        Category = 'Context'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'ReplaceImage'
        HidePluginInDisplay = false
        FlattenPreviewImage = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = true
        Visibility = 'Developer'

        MemoryCachePolicy = 'Permanent'
        DiskCachePolicy = 'Permanent'
    end

    methods (Static)
        function results = RunPlugin(dataset, reporting)
            roi = dataset.GetResult('PTKLungROI');
            left_and_right_lung_mask = dataset.GetResult('PTKLeftAndRightLungs');
            results = PTKGetLeftLungROIFromLeftAndRightLungs(roi, left_and_right_lung_mask, reporting);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function y = apply_3d_func(f,x)

% apply_3d_func - apply a function to 3D arrays
%
%   y = apply_3d_func(f,x);
%
%   f should take a 2D NxNxN image and ouput a 2D NxNxN image
%   x is an N^3 x Q input
%   x is an N^3 x Q input
%
%   This function is used to adapt 3D functions (e.g. filtering) for the
%   barycenter code that assume 1D vectors.
%
%   Copyright (c) 2014 Gabriel Peyre

N = round( (size(x,1))^(1/3) );
P = size(x,2);

resh = @(x)reshape(x,[N N N]);
flat = @(x)x(:);

y = zeros(size(x));
for i=1:P
   y(:,i) = flat(f(resh(x(:,i))));
end

end

author: gpeyre


text:
function vol=smoothbinvol(vol,layer)
%
% vol=smoothbinvol(vol,layer)
%
% perform a memory-limited 3D image smoothing
%
% author: Qianqian Fang <fangq at nmr.mgh.harvard.edu>
%
% input:
%     vol: a 3D volumetric image to be smoothed
%     layer: number of iterations for the smoothing
%
% output:
%     vol: the volumetric image after smoothing
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

dim=size(vol);
dxy=dim(1)*dim(2);
fulllen=prod(dim);

weight=1./6.;
step=4000;

% in case vol is a logical
vol=double(vol);
offs=[1,-1,dim(1), -dim(1),dxy, -dxy];

for i=1:layer
  % find all non-zero values
  idx=find(vol);
  % get the neighbors of all the non-zero values
  % this may cause wrapping -- TODO
  val=vol(idx);
  for k=1:6
    nextidx=idx+offs(k);
    % find all 1-valued voxels that are located within the domain
	  goodidx=find(nextidx>0 & nextidx<fulllen);
    % for all neighboring voxels, add a fraction from the non-0 voxels
    % problematic when running in parallel (racing)
    len=length(goodidx);
    % control granualarity with step
    if(len>step)
        for j=1:step:len-step
            vol(nextidx(goodidx(j:j+step-1)))=vol(nextidx(goodidx(j:j+step-1)))+weight*val(goodidx(j:j+step-1));
        end
        vol(nextidx(goodidx(j+step:end)))=vol(nextidx(goodidx(j+step:end)))+weight*val(goodidx(j+step:end));
    else
        vol(nextidx(goodidx))=vol(nextidx(goodidx))+weight*val(goodidx);
    end
    % the above line may change the values of the non-zero voxels, recover
    % them
  end
  vol(idx)=val;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [attnInfo] = attnLayerForward(h_t, params, model, attnData, maskInfo)
%
% Attentional Layer: from lstm hidden state to softmax hidden state.
% Input:
%   attnData: require attnData.srcHidVecsOrig and attnData.srcLens
%
% Thang Luong @ 2015, <lmthang@stanford.edu>
%
  attnInfo = [];
  if params.attnGlobal % global
    srcHidVecs = attnData.srcHidVecsOrig;
    % attnInfo.srcMaskedIds = [];
    attnInfo.srcMaskedIds = find(attnData.srcMask(:, 1:params.numSrcHidVecs)'==0); % numSrcHidVecs * curBatchSize
  else % local
    [mu, attnInfo] = regressPositions(model, h_t, attnData.srcLens, params);
    srcPositions = floor(mu);

    % assert
    if params.assert
      assert(isempty(find(srcPositions<1,1)));
      assert(isempty(find(srcPositions(maskInfo.unmaskedIds)>(attnData.srcLens(maskInfo.unmaskedIds)-1),1)));
    end

    % reverse
    if params.isReverse
      srcPositions = params.srcMaxLen - srcPositions;
    end

    % build context vectors
    [srcHidVecs, attnInfo] = buildSrcVecs(attnData.srcHidVecsOrig, srcPositions, maskInfo, attnData.srcLens, params.srcMaxLen, params, attnInfo);

    attnInfo.srcMaskedIds = find(attnInfo.alignMask==0);
  end % end else if attnGlobal

  % compute alignScores: numAttnPositions * curBatchSize
  % TODO: precompute for attnOpt2 and attnOpt3 (we can premultiply srcHidVecs with W_a (attnOpt2) or W_a_src (attnOpt3)
  if params.attnOpt==1 || params.attnOpt==2 % dot product or general dot product
    if params.attnOpt==1 % dot product
      [alignScores] = srcCompareLayerForward(srcHidVecs, h_t, params);
    elseif params.attnOpt==2 % general dot product
      attnInfo.transform_ht = model.W_a * h_t; % TODO: shift the multiplication to srcHidVecs
      [alignScores] = srcCompareLayerForward(srcHidVecs, attnInfo.transform_ht, params);
    end
  elseif params.attnOpt==3 % Bengio's style
    % f(H_src + W_a*h_t): lstmSize * (curBatchSize * numAttnPositions))
    attnInfo.src_ht_hid = reshape(params.nonlinear_f(bsxfun(@plus, srcHidVecs, model.W_a*h_t)), params.lstmSize, []);

    % v_a * src_ht_hid
    alignScores = linearLayerForward(model.v_a, attnInfo.src_ht_hid); % 1 * (curBatchSize * numAttnPositions)
    alignScores = reshape(alignScores, params.curBatchSize, params.numAttnPositions)'; % numAttnPositions * curBatchSize
  end

  if params.attnGlobal == 0 && params.normLocalAttn % new approach for local attention after EMNLP'15
    [attnInfo.distWeights, attnInfo.scaleX] = distLayerForward(mu, attnInfo, params); % numAttnPositions*curBatchSize
    attnInfo.unNormAlignWeights =  alignScores .* attnInfo.distWeights; % weighted by distances
    attnInfo.alignWeights = normLayerForward(attnInfo.unNormAlignWeights, attnInfo.srcMaskedIds);
    attnInfo.alignScores = alignScores;
  else
    % normalize -> alignWeights
    attnInfo.alignWeights = normLayerForward(alignScores, attnInfo.srcMaskedIds);

    % local, regression, multiply with distWeights
    if params.attnGlobal == 0
      [attnInfo.distWeights, attnInfo.scaleX] = distLayerForward(mu, attnInfo, params); % numAttnPositions*curBatchSize
      attnInfo.preAlignWeights = attnInfo.alignWeights;
      attnInfo.alignWeights = attnInfo.preAlignWeights.* attnInfo.distWeights; % weighted by distances
    end
  end

  % assert
  if params.assert
    assert(computeSum(attnInfo.alignWeights(attnInfo.srcMaskedIds), params.isGPU)==0);
  end

  attnInfo.alignWeights(:, maskInfo.maskedIds) = 0;
  % alignWeights, srcHidVecs -> contextVecs
  [contextVecs] = contextLayerForward(attnInfo.alignWeights, srcHidVecs, maskInfo.unmaskedIds, params);

  % f(W_h*[context_t; h_t])
  attnInfo.input = [contextVecs; h_t];
  attnInfo.h_t = h_t;
  attnInfo.softmax_h = hiddenLayerForward(model.W_h, attnInfo.input, params.nonlinear_f);

  % assert
  if params.assert
    assert(isequal(size(attnInfo.alignWeights), [params.numAttnPositions, params.curBatchSize]));
    assert(isequal(size(h_t), size(contextVecs))); % lstmSize * curBatchSize
  end
end

function [mu, h2sInfo] = regressPositions(model, h_t, srcLens, params)
  % h_t -> scales=sigmoid(v_pos*f(W_pos*h_t)) in [0, 1]
  [h2sInfo.scales, h2sInfo.posForwData] = scaleLayerForward(model.W_pos, model.v_pos, h_t, params);

  % scales -> srcPositions
  mu = h2sInfo.scales.*(srcLens-1) + 1;
end

author: lmthang


text:
function addHelpSearch

%% Make help files searchable

mfilePath = mfilename('fullpath');
gibbonPath=fileparts(fileparts(mfilePath));
gibbonHelpPath=fullfile(gibbonPath,'docs','html');
builddocsearchdb(gibbonHelpPath);

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function MRSIStruct = op_CSICombineTime(MRSIStruct, dimensionLabel)
    newTimeLength = prod(getSizeFromDimensions(MRSIStruct, {'t', dimensionLabel}));
    [MRSIStruct, prevPermute, prevSize]= reshapeDimensions(MRSIStruct, {'t', dimensionLabel});
    data = getData(MRSIStruct);

    data = reshape(data, newTimeLength, []);
    MRSIStruct = setData(MRSIStruct, data);
    prevSize(2) = [];
    prevSize(1) = newTimeLength;

    prevPermute = removeDimPrevPermute(prevPermute, 2);
    MRSIStruct = reshapeBack(MRSIStruct, prevPermute, prevSize);
end

author: CIC-methods

repo: FID-A

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


text:
% MAIN - Point Mass
%
% Finds the optimal trajectory to slide a point-mass across a 1d
% frictionless plane, using a variety of cost functions.
%
% Simple force-squared cost function  --  This is easy to optimize
%

clc; clear;
addpath ../../

% User-defined dynamics and objective functions
problem.func.dynamics = @(t,x,u)( dynamics(x,u) );
problem.func.pathObj = @(t,x,u)( obj_forceSquared(u) );

% Problem bounds
problem.bounds.initialTime.low = 0;
problem.bounds.initialTime.upp = 0;
problem.bounds.finalTime.low = 1.0;
problem.bounds.finalTime.upp = 1.0;

problem.bounds.state.low = [0; -inf];
problem.bounds.state.upp = [1; inf];
problem.bounds.initialState.low = [0;0];
problem.bounds.initialState.upp = [0;0];
problem.bounds.finalState.low = [1;0];
problem.bounds.finalState.upp = [1;0];

problem.bounds.control.low = -50; %-inf;
problem.bounds.control.upp = 50; %inf;

% Guess at the initial trajectory
problem.guess.time = [0,1];
problem.guess.state = [0, 0; 1, 0];
problem.guess.control = [1, -1];

% Options for fmincon
problem.options.nlpOpt = optimset(...
    'Display','iter',...
    'GradObj','on',...
    'GradConstr','on',...
    'DerivativeCheck','off');   %Fmincon automatically checks

problem.options.method = 'trapezoid';
% problem.options.method = 'rungeKutta';

% Solve the problem
soln = optimTraj(problem);
t = soln.grid.time;
q = soln.grid.state(1,:);
dq = soln.grid.state(2,:);
u = soln.grid.control;

% Plot the solution:
figure(1); clf;

subplot(3,1,1)
plot(t,q)
ylabel('pos')
title('Move Point Mass');

subplot(3,1,2)
plot(t,dq)
ylabel('vel')

subplot(3,1,3)
plot(t,u)
ylabel('force')

author: MatthewPeterKelly

repo: OptimTraj


text:
function used = getvariables(F)

F = flatten(F);
if length(F.clauses) == 0
    used = [];
    return
end

if isa(F.clauses{1},'cell')
    F = flatten(F);
end

used = recursivegetvariables(F,1,length(F.clauses));

function used = recursivegetvariables(F,startindex,endindex)

if endindex-startindex>50
    newstart = startindex;
    mid = ceil((startindex + endindex)/2);
    newend = endindex;
    used1 = recursivegetvariables(F,newstart,mid);
    used2 = recursivegetvariables(F,mid+1,newend);
    used = uniquestripped([used1 used2]);
else
    used = [];
    if startindex <= length(F.clauses)

        if F.clauses{startindex}.type == 56
            used = [];
            for j = 1:length(F.clauses{startindex}.data)
                used = [used getvariables(F.clauses{startindex}.data{j})];
            end
        else
            used = getvariables(F.clauses{startindex}.data);
        end

        for i = startindex+1:endindex

            %Fivars = getvariables(F.clauses{i}.data);
            if F.clauses{i}.type == 56
                % Meta constraint such as implies. This object is just holding
                % the data involved
                Fivars = [];
                for j = 1:length(F.clauses{i}.data)
                    Fivars = [Fivars getvariables(F.clauses{i}.data{j})];
                end
            else
                Fivars = getvariables(F.clauses{i}.data);
            end

            if ~isequal(used,Fivars(:)')
                used = [used Fivars(:)'];
            end
        end
        used = uniquestripped(used);
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [km3] = oz2km3(oz)
% Convert volume from US liquid ounces to cubic kilometers.
% Chad Greene 2012
km3 = oz*2.9573529563e-14;

author: Sable

repo: mcbench-benchmarks


text:
function [Vdraw, phi_Vdraw]=sampleV(Psi,phi_V,Vold,priorValues,dataValues)
% This function takes a draw from the conditional posterior of the variance
% V of the SV in the local mean.

%% Initialize
offset_c=priorValues.offset_c;

p = find(isnan(Psi(:,1)),1,'last');
[T,M] = size(Psi);
% Ppsi=dataValues.Ppsi;

% obtain prior data
startMeanVector=priorValues.mean_ln_v0;
startVarVector=priorValues.var_ln_v0;

priorPhi_V=priorValues.phi_v;
priorD_V  =priorValues.d_v;

%% Prepare Sampling

Vvars=zeros(T,M);
phi_Vdraw=zeros(M,1);

% construct VvarsOld
% VvarsOld=zeros(T,M);
%
% for i=1:M
%     VvarsOld(:,i)=reshape(Vold(i,i,:),T,1,1);
% end

%% Sample Vvars
for i=1:M

    % prepare for each i
        zEntry=p+1;

        PsiDiff=Psi(zEntry+1:T,i)-Psi(zEntry:T-1,i);

        residsTemp=PsiDiff;
        yStar=log(residsTemp.^2+offset_c);
        phi=phi_V(i);
        lnSigma2=log(Vold(zEntry+1:T,i)); % note log => h = ln sigma2

        startMean=startMeanVector(i);
        startVar=startVarVector(i);

        % sample states
        [yStarAdj, Ht] = bear.statesMix(yStar,lnSigma2);

        % sample log variances
        [logVarsDraw_Vi]=bear.KF_CKsimSV(yStarAdj,Ht,phi,startMean,startVar);
        Vvars(:,i)=[0.00000001*ones(zEntry,1) ;exp(logVarsDraw_Vi)];

        % sample phi
        [phiDraw_Vi]=bear.samplePhi(logVarsDraw_Vi,priorD_V,priorPhi_V);
        phi_Vdraw(i)=phiDraw_Vi;

end

%% Construct V
% Vdraw=zeros(M,M,T);
%
% for i=1:M
%     Vdraw(i,:)=Vvars(:,i);
% end
Vdraw=Vvars;
end

author: european-central-bank

repo: BEAR-toolbox


text:
function obj = updateContourProjection(obj,contourIndex)

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

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

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

%-PLOT DATA STRUCTURE- %
contour_data = get(obj.State.Plot(contourIndex).Handle);

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

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

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

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

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

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

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

%-contour name-%
obj.data{contourIndex}.name = contour_data.DisplayName;

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

%-setting the plot-%
xdata = contour_data.XData;
ydata = contour_data.YData;
zdata = contour_data.ZData;

%-contour type-%
obj.data{contourIndex}.type = 'surface';

%-contour x and y data
obj.data{contourIndex}.x = xdata;
obj.data{contourIndex}.y = ydata;

%-contour z data-%
obj.data{contourIndex}.z = zdata;%-2*ones(size(zdata));

%-setting for contour lines z-direction-%
obj.data{contourIndex}.contours.z.start = contour_data.LevelList(1);
obj.data{contourIndex}.contours.z.end = contour_data.LevelList(end);
obj.data{contourIndex}.contours.z.size = contour_data.LevelStep;
obj.data{contourIndex}.contours.z.show = true;
obj.data{contourIndex}.contours.z.usecolormap = true;
obj.data{contourIndex}.hidesurface = true;
obj.data{contourIndex}.surfacecolor = zdata;

obj.data{contourIndex}.contours.z.project.x = true;
obj.data{contourIndex}.contours.z.project.y = true;
obj.data{contourIndex}.contours.z.project.z = true;

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

%-contour visible-%

obj.data{contourIndex}.visible = strcmp(contour_data.Visible,'on');

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

%-contour showscale-%
obj.data{contourIndex}.showscale = false;

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

%-colorscale (ASSUMES PATCH CDATAMAP IS 'SCALED')-%
colormap = figure_data.Colormap;

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

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

%-contour reverse scale-%
obj.data{contourIndex}.reversescale = false;

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

%-aspect ratio-%
ar = obj.PlotOptions.AspectRatio;

if ~isempty(ar)
    if ischar(ar)
        obj.layout.scene.aspectmode = ar;
    elseif isvector(ar) && length(ar) == 3
        xar = ar(1);
        yar = ar(2);
        zar = ar(3);
    end
else

    %-define as default-%
    xar = max(xdata(:));
    yar = max(ydata(:));
    xyar = max([xar, yar]);
    zar = 0.6*xyar;
end

obj.layout.scene.aspectratio.x = xyar;
obj.layout.scene.aspectratio.y = xyar;
obj.layout.scene.aspectratio.z = zar;

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

%-camera eye-%
ey = obj.PlotOptions.CameraEye;

if ~isempty(ey)
    if isvector(ey) && length(ey) == 3
        obj.layout.scene.camera.eye.x = ey(1);
        obj.layout.scene.camera.eye.y = ey(2);
        obj.layout.scene.camera.eye.z = ey(3);
    end
else

    %-define as default-%
    xey = - xyar; if xey>0 xfac = -0.2; else xfac = 0.2; end
    yey = - xyar; if yey>0 yfac = -0.2; else yfac = 0.2; end
    if zar>0 zfac = 0.2; else zfac = -0.2; end

    obj.layout.scene.camera.eye.x = xey + xfac*xey;
    obj.layout.scene.camera.eye.y = yey + yfac*yey;
    obj.layout.scene.camera.eye.z = zar + zfac*zar;
end

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

%-contour showlegend-%

leg = get(contour_data.Annotation);
legInfo = get(leg.LegendInformation);

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

obj.data{contourIndex}.showlegend = showleg;

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

end

author: plotly


text:
function joined = strjoin(tokens, delimiter)
% strjoin Joins multiple strings
%
% Joins a cell array of strings in a single string with a given delimiter in between.
%
% Input:
% - tokens (cell): A cell array of strings.
% - delimiter (string): A delimiter string.
%
% Output:
% - joined (string): Joined string.
%

if nargin < 2
    delimiter = ' ';
end;

if ~iscell(tokens)
    joined = tokens;
    return;
end;

if isempty(tokens)
    joined = '';
    return;
end;

tokens = tokens(:)';

if numel(tokens) == 1
    joined = tokens{1};
    return;
end;

joined = cellfun(@(x) [delimiter, x], tokens(2:end), 'UniformOutput', false);
joined = [tokens{1} joined{:}];

author: votchallenge

repo: toolkit-legacy


text:
function [mu, Q] = matrixCoherence(A)
% Computes the coherence of a matrix `A`, that is
% the maximum absolute value of the cross correlations
% between the columns of `A`
%
% USAGE:
%
%    [mu, Q] = matrixCoherence(A)
%
% INPUT:
%    A:     Matrix
%
% OUTPUTS:
%    mu:    coherence
%    Q:     matrix used by `mu`

n = size(A, 2);
Q = zeros(n, n);
for j = 1:n
    for k = 1:n
        if j ~= k
            Q(j, k) = abs(A(:, j)' * A(:, k)) / (norm(A(:, j)) * norm(A(:, k)));
        end
    end
end
mu = max(max(Q));

author: opencobra

repo: cobratoolbox


text:
function Q=quatMult(quat2,quat1,handed)
%%QUATMULT  Multiply two quaternions together. Quaternions are an extension
%           of complex numbers. The type of product performed here is
%           deemed the Grassman product. Unit-magnitude quaternions are
%           often used to represent rotations and the multiplication
%           quat2*quat1 is the orientation obtained by rotating by quat1
%           and then rotating by quat2. Quaternion multiplication is not
%           commutative. The quaternions can be either left or
%           right-handed.
%
%INPUTS:quat2 A 4XN set of N quaternions, where the first element in each
%             column is the scalar part of the quaternion (sometimes called
%             q0 or q4) and the next three elements are the (hypercomplex)
%             vector part. That is, the hypercomplex quaternion given by
%             quat2(:,1) can be written in hypercomplex, non-vector form
%             as quat2(1,1)+i*quat2(2,1)+j*quat2(3,1)+k*quat2(4,1), where
%             i, j, and k are all roots of -1. When considering unit
%             quaternions, this means that it has the form
%             [cos(theta/2);sin(theta/2)u'], where theta is an angle and u
%             is a 3X1 unit vector.
%       quat1 A 4XN set of N quaternions that are to be left-multiplied by
%             the corresponding quaternions in quat2.
%     handed The handedness of the quaternions. If omitted or an empty
%            matrix is passed,, it is assumed that the quaternions are
%            right-handed (the standard). Possible values are
%            'right' The default if omitted. The quaternion multiplication
%                    is assumed right-handed (standard).
%            'left'  The quaternion multiplication is assumed left-handed.
%                    This is used in someplaces, including the reference
%                    from Shuster, below.
%
%OUTPUTS: Q   The 4XN matrix of quaternions products such that
%             Q(:,i)=quat2(:,i)*quat1(:,1), where the multiplication
%             operation is quaternion multiplication.
%
%Properties of quaternions including multiplication are described in [1].
%Details of rotation using unit quaternions are given in [2].
%
%The difference in multiplication between right-handed and left-handed
%quaternions lies in the sign used for the cross product.
%
%A quaternion of the form q(1)+i*q(2)+j*q(3)+k*q(4) that obeys right-handed
%multiplication rules supports the following rules for multiplication of i,
%j, and k, where an item in a row is multiplied by an item in the column to
%get the result:
%  i,  j, k
%i -1, k,-j
%j -k,-1, i
%k  j,-i,-1
%On the other hand, left-handed multiplication rules flip the signs of the
%off-diagonal terms:
%  i,  j, k
%i -1,-k, j
%j  k,-1,-i
%k -j, i,-1
%
%REFERENCES:
%[1] Weisstein, Eric W. "Quaternion." From MathWorld--A Wolfram Web
%    Resource. http://mathworld.wolfram.com/Quaternion.html
%[2] M. D. Shuster, "A survey of attitude representations," The Journal of
%    the Astronautical Sciences, vol. 41, no. 4, pp. 439-517, Oct.-Dec.
%    1993.
%
%September 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%The number of quaternions.
N=size(quat2,2);

if(nargin<3||isempty(handed))
    handed='right';
end

Q=zeros(4,N);
%The handedness switches the sign of the cross product
switch(handed)
    case 'right'
        for curQuat=1:N
            s1=quat1(1,curQuat);
            v1=quat1(2:4,curQuat);
            s2=quat2(1,curQuat);
            v2=quat2(2:4,curQuat);

            Q(:,curQuat)=[s1*s2-dot(v1,v2);
                          s1*v2+s2*v1+cross(v2,v1)];
        end
    case 'left'
        %Flip the algebra to left-handed for the formula below.
         for curQuat=1:N
            s1=quat1(1,curQuat);
            v1=quat1(2:4,curQuat);
            s2=quat2(1,curQuat);
            v2=quat2(2:4,curQuat);

            Q(:,curQuat)=[s1*s2-dot(v1,v2);
                          s1*v2+s2*v1-cross(v2,v1)];
        end
    otherwise
        error('Invalid handedness provided.')
end
end

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

author: USNavalResearchLaboratory


text:
function updateContents(folder)
%UPDATECONTENTS Create a Contents.m file including subdirectories
%
%   IOSR.GENERAL.UPDATECONTENTS scans through the current directory, and
%   its subdirectories, and builds a Contents file similar to Matlab's
%   report-generated Contents.m files. Any existing Contents.m file will be
%   overwritten.
%
%   IOSR.GENERAL.UPDATECONTENTS(FOLDER) scans through the directory FOLDER.
%
%   Typing
%
%       help(FOLDER)
%
%   or
%
%       help path/to/folder
%
%   will display Contents.m in the Command Window, and display links to the
%   help for any functions that are in Matlab's search path.
%
%   NB: Do not use Matlab's Contents Report generator to edit the
%   Contents.m file. Execute this function to update it.

%   Copyright 2016 University of Surrey.

    % apply function in current directory
    if nargin<1
        folder = cd;
    end

    % check input is valid
    assert(ischar(folder), 'iosr:updateContents:invalidPath', '''directory'' should be a charater array (string)')
    assert(exist(folder,'dir')==7, 'iosr:updateContents:invalidPath', [folder ' does not exist'])

    % check last character in path is not filesep (e.g. '/')
    if strcmp(folder(end),filesep)
        folder = folder(1:end-1);
    end
    fIX = strfind(folder,filesep); fIX = fIX(end);

    % name of file to create
    filename = 'Contents.m';

    % Name of the folder
    [~,name] = fileparts(folder);

    % delete if it already exists
    if exist([folder filesep filename],'file')==2
        delete([folder filesep filename])
    end

    % get subfolders
    dirs = iosr.general.getContents(folder,'filter','folders','rec',true,'path','full','sort',true);
    dirs = [{folder}; dirs];

    % get files
    files = cell(0,1);
    H1_lines = cell(0,1);
    for d = 1:length(dirs)
        temp = iosr.general.getContents(dirs{d},'filter','files','sort',true);
        if ~isempty(temp)
            temp = temp(cellfun(@(x) isempty(strfind(x,'~')),temp)); % remove temporary files
            temp = temp(cellfun(@(x) isempty(strfind(x,'.mex')),temp)); % remove compiled mex files
            temp = temp(cellfun(@(x) ~strcmp(x,filename),temp)); % remove Contents.m
            H1_lines = [H1_lines; {''}; {''}]; %#ok<AGROW> % insert blank lines where no functions will be
            % determine package prefix
            pkgprefix = strrep(dirs{d},[filesep '+'],'.');
            pkgprefix = strrep(pkgprefix,[filesep '@'],'.');
            dots = strfind(pkgprefix,'.');
            if ~isempty(dots)
                pkgprefix = [pkgprefix(dots(1)+1:end) '.'];
            else
                pkgprefix = '';
            end
            for f = 1:length(temp) % read H1 lines
                H1_lines = [H1_lines; {get_H1_line([dirs{d} filesep temp{f}])}]; %#ok<AGROW> % add H1 lines
                % remove extension from and add package prefix to m-files
                [~,fname,ext] = fileparts(temp{f});
                if strcmpi(ext,'.m')
                    temp{f} = [pkgprefix fname];
                end
            end
            files = [files; {''}; {upper(dirs{d}(fIX+1:end))}; temp;]; %#ok<AGROW> % add filenames
        end
    end

    % longest file name (so appropriate space can be added between files and H1 lines
    longest_word = max(cellfun(@length,files(cellfun(@(x) ~isempty(x),H1_lines))));

    % write to output
    nrows = length(files);
    fid = fopen(filename, 'w'); % open file for writing
    fprintf(fid, '%s\n%% \n', ['% ' upper(name)]);
    fprintf(fid, '%s\n', ['%   Contents file for ' upper(folder(fIX+1:end)) ' and its subfolders.']);
    for row=1:nrows
        if isempty(H1_lines{row})
            fprintf(fid, '%s\n', ['%   ' files{row,:}]);
        else
            rowfilename = files{row,:};
            [~,name,ext] = fileparts(rowfilename);
            if strcmpi(ext,'.m') % remove extension from m files
                rowfilename = name;
            end
            fprintf(fid, '%s\n',['%   ' rowfilename repmat(' ',1,longest_word-length(rowfilename)) ' - ' H1_lines{row,:}]);
        end
    end
    fprintf(fid, '%%    \n%%   %s on %s at %s.\n', 'This file was generated by updateContents.m',datestr(datetime('now'),'dd mmm yyyy'),datestr(datetime('now'),'HH:MM:SS'));
    fclose(fid);

end

function H1_line = get_H1_line(filename)
%GET_H1_LINE get the H1 line for a file

    [~,name,ext] = fileparts(filename);
    H1_line = ''; % default output
    if strcmp(ext,'.m')
        fid = fopen(filename); % open file
        tline = fgetl(fid); % read first line
        while ischar(tline)
            k = strfind(tline,'%'); % find comment
            if ~isempty(k) % if it is found
                k = k(1);
                ispercents = false(size(tline(k:end)));
                ispercents(strfind(tline(k:end),'%'))=true;
                start = k+find(~(isspace(tline(k:end)) | ispercents),1,'first')-1;
                if ~isempty(start)
                    tline = tline(start:end); % remove leading space/percent
                    IX = strfind(lower(tline),lower(name));
                    if ~isempty(IX)
                        if IX(1)==1
                            tline = tline(length(name)+1:end); % remove function name
                        end
                        tline = strtrim(tline); % remove any leading/trailing space
                    end
                    H1_line = tline;
                    H1_line = strtrim(H1_line);
                    if ~isempty(H1_line)
                        if strcmp(H1_line(end),'.') % remove trailing period
                            H1_line = H1_line(1:end-1);
                        end
                        H1_line(1) = upper(H1_line(1)); % capitalize first letter
                    end
                end
                tline = -1; % set tline to numeric
            else
                tline = fgetl(fid); % read next line
            end
        end
        fclose(fid);
    end

end

author: IoSR-Surrey

repo: MatlabToolbox


text:
function HDF5_DatasetWriteStrings(fname, location, data)

if ~exist(fname, 'file')
    fid = H5F.create(fname, 'H5F_ACC_TRUNC', 'H5P_DEFAULT', 'H5P_DEFAULT');
else
    fid = H5F.open(fname,'H5F_ACC_RDWR','H5P_DEFAULT');
end
data = cell2str_new(data);

filetype = H5T.copy('H5T_FORTRAN_S1');
H5T.set_size(filetype, size(data,2));
memtype = H5T.copy('H5T_C_S1');
H5T.set_size(memtype, size(data,2));

% Create dataspace. Setting maximum size to [] sets the maximum
% size to be the current size.
space = H5S.create_simple(1, size(data,1), []);

% Create the dataset and write the string data to it.
dset = H5D.create(fid, location, filetype, space, 'H5P_DEFAULT');

% Transpose the data to match the layout in the H5 file to match C
% generated H5 file.
H5D.write(dset, memtype, 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT',  data');

% Close and release resources.
H5D.close(dset);
H5S.close(space);
H5T.close(filetype);
H5T.close(memtype);
H5F.close(fid);


text:
function [lat, long ] = gcj2wgs_exact(gcjLat, gcjLng)
% adapted from https://github.com/googollee/eviltransform

    initDelta = 0.01;
    threshold = 0.000001;
    dLat = initDelta;
    dLng = initDelta;
    mLat = gcjLat - dLat;
    mLng = gcjLng - dLng;
    pLat = gcjLat + dLat;
    pLng = gcjLng + dLng;

    for i=1:30
        wgsLat = (mLat + pLat) / 2;
        wgsLng = (mLng + pLng) / 2;
        [tmplat, tmplng] = wgs2gcj(wgsLat, wgsLng);
        dLat = tmplat - gcjLat;
        dLng = tmplng - gcjLng;
        if all(abs(dLat) < threshold & abs(dLng) < threshold)
            lat = wgsLat;
            long =wgsLng;
            return;
        end
        pLat(dLat>=0) = wgsLat(dLat>=0);
        mLat(dLat<=0) = wgsLat(dLat<=0);
        pLng(dLng>=0) = wgsLng(dLng>=0);
        mLng(dLng<=0) = wgsLng(dLng<=0);
    end

   lat = wgsLat;
   long =wgsLng;

end

author: googollee

repo: eviltransform


text:
function [hklGrid,wGrid,id] = gridify(hkl,varargin)
% approximate a list of Miller indices with Miller indices on a grid
%
% Syntax
%   [hklGrid, weightsG, id] = gridify(vec)
%   [hklGrid, weightsG, id] = gridify(vec,'weights',weights)
%
% Input
%  vec - @vector3d
%
% Output
%  hklGrid - @Miller
%  weightsG - double
%  id    -
%
% Options
%  weights    -
%  resolution -

sR = hkl.CS.fundamentalSector;
hkl = project2FundamentalRegion(hkl);

[hklGrid,wGrid,id] = gridify@vector3d(hkl,sR,varargin{:});

hklGrid = Miller(hklGrid,hkl.CS);

author: mtex-toolbox

repo: mtex

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


text:
function loglik = compute_log_lik(P, S_bar, V, E_z, E_zz, RO, c, Tr, sigma_sq)

[K, T] = size(E_z);
J = size(S_bar, 2);

M_t = zeros(2*J, K);

loglik = - 0.5*T * (2*J*log(sigma_sq));
for t = 1:T,
   %%%%%%%%%%%%%%%%%%%%%%%%% Changed code here %%%%%%%%%%%%%%%%%%%%%%
   %R_t = RO{t};
   R_t  = c(t,1)*RO{t};
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   Sdef = S_bar;
   for kk = 1:K,
      Sdef = Sdef + E_z(kk,t)*V((kk-1)*3+[1:3],:);

      M_t(1:J, kk) = (R_t(1,:)*V((kk-1)*3+[1:3],:))';
      M_t(J+1:end, kk) = (R_t(2,:)*V((kk-1)*3+[1:3], :))';
   end;

   invSigmaSq_p = eye(2*J)./sigma_sq;

   f_bar_t = R_t(1:2,:)*S_bar;
   f_bar_t = [f_bar_t(1,:) f_bar_t(2,:)]';

   f_t = [P(t, :) P(t+T, :)]';
   t_vect_t = [Tr(t,1)*ones(J,1); Tr(t,2)*ones(J,1)];

   covZ_t = E_zz((t-1)*K+1:t*K,:) - E_z(:,t)*E_z(:,t)';
   loglik = loglik - 0.5*(((f_t-f_bar_t-t_vect_t)./sigma_sq)'*(f_t-f_bar_t-t_vect_t)) + (((f_t-f_bar_t-t_vect_t)'./sigma_sq)*M_t*E_z(:,t)) ...
      - 0.5*trace(((M_t./sigma_sq)'*M_t) * E_zz((t-1)*K+1:t*K,:)) - 0.5*log(det(covZ_t));
end

repo: CategoryShapes


text:
% testing on new December 2016 price data

dataArray = csvread('results.csv');
prices = transpose(dataArray(:,2));
askVolume = dataArray(:,3);
bidVolume = dataArray(:,4);

prices = prices(1:2:end);
askVolume = askVolume(1:2:end);
bidVolume = bidVolume(1:2:end);

% estimate transaction fee at %1
[error,jinzhi,bank,buy,sell,proba] = brtrade(prices,bidVolume,askVolume, 1);

% set up plots
make_plots(prices, buy, sell, proba, bank, error);

author: panditanvita

repo: BTCpredictor


text:
function c = inf_(a)
%INF_         Implements  inf(a)  (cures problems with inf)
%
%   c = inf(a)
%
% On return, inf(a) <= alpha for all entries alpha in a
%
%************************************************************************
%********  due to conflict with internal variable inf (infinity)  *******
%********                    use function inf_                    *******
%************************************************************************
%

% written  10/05/99     S.M. Rump
% modified 09/02/00     S.M. Rump  rounding unchanged after use
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    take care of huge arrays
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/10/07     S.M. Rump  performance, huge arrays
%

  if a.complex
    if isequal(a.rad,0)                  % faster for sparse matrices
      c = a.mid;
    else
      e = 1e-30;
      if 1+e==1-e                        % fast check for rounding to nearest
        rndold = 0;
      else
        rndold = getround;
      end
      setround(-1)
      if issparse(a.rad)
        [m,n] = size(a.rad);
        [I,J,arad] = find(a.rad);
        c = a.mid - sparse(I,J,complex(arad,arad),m,n);
      else
        c = a.mid - complex(a.rad,a.rad);
      end
      setround(rndold)
    end
  else
    c = a.inf;
  end

repo: OpenMAS


text:
function name = getComputerName()
% GETCOMPUTERNAME returns the name of the computer (hostname)
% name = getComputerName()
%
% WARN: output string is converted to lower case
%
%
% See also SYSTEM, GETENV, ISPC, ISUNIX
%
% m j m a r i n j (AT) y a h o o (DOT) e s
% (c) MJMJ/2007
% MOD: MJMJ/2013

[ret, name] = system('hostname');

if ret ~= 0,
   if ispc
      name = getenv('COMPUTERNAME');
   else
      name = getenv('HOSTNAME');
   end
end
% name = strtrim(lower(name));

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
 function ys = ifftn_fast(xs)
%function ys = ifftn_fast(xs)
%|
%| For some reason, matlab's ifftn routine is suboptimal
%| for the case of 2D FFTs, at least on some machines.
%| The improvement herein was found by Hugo Shi.
%|
%| Note: matlab's ifft() and ifftn() handle an optional second "N" argument in
%| different ways!  So to be safe I am not allowing any second argument here.
%|
%| Copyright 2004-6-28, Jeff Fessler, University of Michigan

if nargin ~= 1, ir_usage, end
if streq(xs, 'test'), ifftn_fast_test, return, end

% around version 7.4, ifftn was fastest, so hardwire that!
ys = ifftn(xs);
return

if ndims(xs) == 2 % 2D or 1D cases
	if min(size(xs)) == 1 % 1D
		ys = ifft(xs);
	else
		ys = ifftn_fast_fftfft(xs);
	end
else
	ys = ifftn(xs);
end

function ys = ifftn_fast_fftfft(xs)
ys = ifft(ifft(xs).').';

% test configuration of ifftn_fast for this machine
function ifftn_fast_test
ifftn_fast_test2
%ifftn_fast_test3

% test configuration of ifftn_fast for this machine for 2D
function ifftn_fast_test2

n = 2^9;
rng(0)
x = rand(n,n);

printm('starting test; be patient.')

% first loop is to get everything in cache or whatever.
% doing it twice is the only way to get an accurate comparison!
for nloop = [2 40];
	tic, for ii=1:nloop, y{1} = ifftn_fast(x); end
	tt(1) = toc; ty{1} = 'fftn_fast';

	tic, for ii=1:nloop, y{2} = ifftn(x); end
	tt(2) = toc; ty{2} = 'fftn';

	tic, for ii=1:nloop, y{3} = ifft(ifft(x).').'; end
	tt(3) = toc; ty{3} = 'fftfft_transpose';

	tic, for ii=1:nloop, y{4} = ifft(ifft(x, [], 1), [], 2); end
	tt(4) = toc; ty{4} = 'fftfft_brack';

	tic, for ii=1:nloop, y{5} = ifft2(x); end
	tt(5) = toc; ty{5} = 'ifft2';

	tic, for ii=1:nloop, y{6} = ifftn_fast_fftfft(x); end
	tt(6) = toc; ty{6} = 'fftfft_func';
end

for ii = 1:length(tt)
	printm('time %19s = %g', ty{ii}, tt(ii))
	if max_percent_diff(y{1}, y{ii}) > 1e-11, error 'bug', end
end

printm('ifftn / ifftn_fast = %g%% ', tt(2) / tt(1) * 100.)
if tt(1) > 1.40 * min(tt(2:end))
	warn 'ifftn_fast is configured supoptimally for your machine!'
else
	printm('ifftn_fast is configured appropriately for your machine')
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function meanRot=meanRotation(param1,w)
%%MEANROTATION Perform a weighted or unweighted average of a set of
%              rotation matrices or unit quaternions (which represent
%              directions). The weights are scalar and the cost function
%              finds the rotation matrix REst to minimize the sum of
%              w(i)*norm(REst-R{i},'fro'), where R{i} is the ith rotation
%              matrix to be averaged (either given explicitly as a
%              rotation matrix or given as a quaternion).
%
%INPUTS: param1 A 3X3XN hypermatrix of N rotation matrices that are to be
%               averaged or a 4XN set of unit-quaternions that are to be
%               averaged. The function automatically detects whether
%               rotation matrices or quaternions were passed. The
%               handedness of the quaternions does not matter.
%             w An optional set of N weights for averaging the rotations.
%               If omitted, all of the rotation matrices/ unit quaternions
%               in param1 are equally weighted.
%
%OUTPUTS: meanRot The weighted mean 3X3 rotation matrix or 4X1 rotation
%                 quaternion, depending on the data type passed in param1.
%
%The algorithm of [1] for performing a weighted average of quaternions or
%rotation matrices is used. Note that the paper contains an additional
%technique for averaging quaternions with non-scalar weights, which is not
%implemented here.
%
%A quaternion can represent an orientation (rotation) in space as a 4X1
%unit vector. However, one cannot simply take the weighted average of a
%set of quaternions to get an average quaternion (that is then normalized)
%to represent the average rotation, because the quaternion representation
%is not unique. Thus, while q and -q represent the same rotation, they
%affect such an average differently, which is bad. This is why the more
%complicated algorithms of the above paper are used.
%
%REFERENCES:
%[1] F. L. Markley, Y. Cheng, J. L. Crassidis, and Y. Oshman, "Averaging
%    quaternions," Journal of Guidance, Control, and Dynamics, vol. 30,
%    no. 4, pp. 1193-1196, Jul. - Aug. 2007.
%
%September 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(size(param1,1)==3)
    usingQuaternions=false;
    RotMats=param1;

    numVals=size(param1,3);
else
    usingQuaternions=true;
    q=param1;

    numVals=size(param1,2);
end

if(nargin<2)
   w=ones(numVals,1);
end

%If quaternions are being averaged rather than rotation matrices.
if(usingQuaternions)
    %Build M as in Equation 2 of [1].
    M=zeros(4,4);
    for curQuat=1:numVals
        M=M+w(curQuat)*(q(:,curQuat)*q(:,curQuat)');
    end

    %The solution to the cost function in Equation 13 is the eigenvector
    %corresponding to the largest eigenvalue in M.
    [V,D]=eig(M);
    %Get the indices of the sorted eigenvalues (and thus of the
    %corresponding eigenvectors).
    [~,idx]=sort(diag(D),'descend');
    %The rotation quaternion is the eigenvector associated with the
    %largest eigenvalue.
    meanRot=V(:,idx(1));%A quaternion
else%If rotation matrices are being averaged rather than quaternions.
    %Build B as in Equation 5 in [1].
    B=zeros(3,3);
    for curRot=1:numVals
       B=B+w(curRot)*RotMats(:,:,curRot);
    end

    %This solves the equivalent formulation of Wahba's problem for
    %Equation 4.
    meanRot=findTransParam(B,eye(3));%A rotation matrix.
end

end

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

repo: TrackerComponentLibrary


text:
% Matlab script for RAPID equivalent commands
% please note that variables should be
% declared as global as in RAPID syntaxis
% this script demonstrates the use of functions in the matlab script. This
% script can be translated to RAPID by means of the matlab2RAPID script.
%
% IN ORDER TO SIMULATE THE PROGRAM:
%   A) FIRST, LOAD A ROBOT
%       robot = load_robot('YASKAWA','yaskawa_motoman_es165d_100');
%
%   B) NEXT, LOAD SOME EQUIPMENT.
%       robot.equipment{1} = load_robot('equipment','conveyor_belt1');
%
%   C) NOW, LOAD AN END TOOL
%       robot.tool= load_robot('equipment','end_tools/vacuum_1');
%   D) FINALLY, LOAD A PIECE TO GRAB BY THE ROBOT
%       robot.piece=load_robot('equipment','aluminum_plate1');
%
%   E) IF NECESSARY, CHANGE THE POSITION AND ORIENTATION OF THE piece,
%   relative to the robot's base reference system.
%
%       robot.piece.T0= [1 0 0 -0.1;
%                        0 1 0 -0.5;
%                        0 0 1 0.2;
%                        0 0 0 1];
%
% during the simulation, call simulation_open_tool; to open the tool and
% simulation_close_tool; to close it.
%
% To grip the piece, call simulation_grip_piece; and
% simulation_release_piece to release it.
%
% The call to each function must be in a correct, otherwise, the result may be
% undesired, thus, typically the correct sequence is:
% simulation_open_tool;
% (approach the piece to grab with MoveL or MoveJ)
% simulation_close_tool;
% simulation_grip_piece; --> the piece will be drawn with the robot
%
% move to a different place with MoveJ or MoveL
% simulation_open_tool;
% simulation_release_piece
%
% In RAPID, consider the use of ConfJ on, or ConfJ off, or ConfL on, or ConfL off
% in the case the controller avoids the change of configurations between target points

function robot_test
global robot TD_vacuum RT_initial RT_approach1 RT_vacuum_catch RT_approach2 RT_vacuum_release RT_transition1 RT_transition2 RT_transition3

%Comment the following lines to avoid loading the robot at every simulation
robot = load_robot('MOTOMAN','ES165D_100');
robot.equipment{1} = load_robot('equipment','conveyor_belt1');
robot.tool= load_robot('equipment','end_tools/vacuum_1');
robot.piece{1}=load_robot('equipment','aluminum_plate1');

adjust_view(robot);
fprintf('Press any key to continue')
%robot.tool.piece_gripped=0;
drawrobot3d(robot, robot.q);

%define the tool
%In RAPID this is done by means of the tooldata structure
TD_vacuum=[1,[[0,0,0.4],[1,0,0,0]],[0.1,[0,0,0.100],[1,0,0,0],0,0,0]];

%define target points FOR SIMULATION
RT_initial=[[2.5250, -0.0000, 2.0500],[0.7071, -0.0000, 0.7071, -0.0000], [0, -1, 2, 2], [9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
RT_approach1=[[2.2000, 0.0000, 1.0000],[0.0000, 0.7071, 0.7071, -0.0000], [0, -1, 1, 0], [9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
RT_vacuum_catch=[[2.2000, 0.0000, 0.8510],[0.0000, 0.7071, 0.7071, -0.0000], [0, -1, 1, 0], [9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
RT_transition1=[[2.2000, 0.0000, 1.3000],[0.0000, -0.0000, 1.0000, -0.0000], [0, -1, 0, 0], [9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
RT_transition2=[[1.1000, -1.4600, 1.3000],[0.0000, 0.0000, 1.0000, 0.0000], [-1, -1, -1, 0], [9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
RT_transition3=[[0.7000, -1.4000, 1.5000],[0.0000, -0.0000, 1.0000, -0.0000], [-1, -1, -1, 0], [9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
RT_approach2=[[0.0000, -1.5000, 0.4000],[0.0000, 0.7071, 0.7071, 0.0000], [-2, -1, 0, 0], [9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];
RT_vacuum_release=[[0.0000, -1.5000, 0.0510],[0.0000, 0.7071, 0.7071, 0.0000], [-2, -1, 0, 0], [9E+09,9E+09,9E+09,9E+09,9E+09,9E+09]];

main;
end

function main()

global TD_vacuum RT_initial RT_approach1 RT_vacuum_catch RT_approach2 RT_vacuum_release RT_transition1 RT_transition2 RT_transition3

%close the tool
simulation_close_tool; %Set do1;

%move to the initial point
MoveJ(RT_initial, 'vmax' , 'fine' ,TD_vacuum, 'wobj0');

% Now open the tool
simulation_open_tool; %Reset do1;

%Move to the approaching point
MoveJ(RT_approach1, 'vmax' , 'fine' ,TD_vacuum, 'wobj0');
%Now, go down to the grabbing target point and
MoveL(RT_vacuum_catch, 'vmax' , 'fine' ,TD_vacuum, 'wobj0');

%and close the tool and grip the piece. These two functions
% must be called to simulate that the gripper has the piece grabbed and the
%tool is closed
simulation_close_tool; %Set do1;
simulation_grip_piece;

%Move to a transition point to secure the integrity of the conveyor belt
MoveL(RT_transition1,'vmax','fine',TD_vacuum,'wobj0');

MoveL(RT_transition2,'vmax','fine',TD_vacuum,'wobj0');

%Move to the approach next to the packaging area
MoveJ(RT_approach2, 'vmax' , 'fine' ,TD_vacuum, 'wobj0');
%go down to the release point inside the box
MoveL(RT_vacuum_release, 'vmax' , 'fine' ,TD_vacuum, 'wobj0');

%yes, release the piece
simulation_open_tool; %Reset do1;
simulation_release_piece; %Reset do1;

%now, go up
MoveL(RT_approach2, 'vmax' , 'fine' ,TD_vacuum, 'wobj0');

%Move to a transition point to secure the integrity of the conveyor belt
MoveJ(RT_transition3,'vmax','fine',TD_vacuum,'wobj0');

%Now, go back to the initial point
MoveJ(RT_initial, 'vmax' , 'fine' ,TD_vacuum, 'wobj0');

end

repo: ARTE


text:
function out = isinf(f)
%ISINF   Test if a CHEBTECH is unbounded.
%   ISINF(F) returns TRUE if F has any infinite values and FALSE otherwise.
%
% See also ISFINITE, ISNAN.

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

% Check if any coefficients are infinite:
out = any(isinf(f.coeffs(:)));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function t1 = joinTrees(t1, t2, alpha)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% t12 = joinTrees(t1, t2, alpha)
%   create a new KD-tree with t1 and t2 as the children of the root
%   The t1 subtree recieves weight alpha; the t2 subtree has wt. 1-alpha
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2003 Alexander Ihler; distributable under GPL -- see README.txt

if(t1.D ~= t2.D)
  error('Input trees have different dimensionality!');
end

if(isfield(t1, 'type') ~= isfield(t2, 'type'))
  error('Cant merge BallTree with BallTreeDensity!');
end

if(t1.type ~= t2.type)
  error('Trees must have the same type of kernel!');
end

if(nargin < 2 || nargin > 3)
  error('wrong number of arguments');
end
if(nargin == 2)
  alpha = 0.5;
end

% leave off one of the zeros between leaves and nodes
if(t1.N==1) t1nodes = []; t1leaves=[1];
  else t1nodes = [1:t1.N-1]; t1leaves = [t1.N+1:2*t1.N]; end;
if(t2.N==1) t2nodes = []; t2leaves = [1];
  else t2nodes = [1:t2.N-1]; t2leaves = [t2.N+1:2*t2.N]; end;

if(t1.N~=1) t1root=2;
else t1root=3+length(t2nodes);
end
if(t2.N~=1) t2root=2+length(t1nodes);
else t2root=3+length(t1nodes)+length(t1leaves);
end
t1N = t1.N;
t1.N = t1.N + t2.N;

Os = zeros(t1.D, 1);

t1.centers = [Os t1.centers(:,t1nodes) t2.centers(:,t2nodes) Os ...
              t1.centers(:,t1leaves) t2.centers(:,t2leaves)];
t1.ranges = [Os t1.ranges(:,t1nodes) t2.ranges(:,t2nodes) Os t1.ranges(:,t1leaves) ...
             t2.ranges(:,t2leaves)];
t1.weights = [0 alpha*t1.weights(t1nodes) (1-alpha)*t2.weights(t2nodes) 0 ...
              alpha*t1.weights(t1leaves) (1-alpha)*t2.weights(t2leaves)];

t1.means = [Os t1.means(:,t1nodes) t2.means(:,t2nodes) Os ...
            t1.means(:,t1leaves) t2.means(:,t2leaves)];

% take care of variable BWs
t1varBWs = size(t1.bandwidth, 2) > 2*t1N;
t2varBWs = size(t2.bandwidth, 2) > 2*t2.N;
varBWs = t1varBWs || t2varBWs || ...
         sum(t1.bandwidth(:,1) ~= t2.bandwidth(:,1)) > 0;
if(varBWs)
  t1.bandwidth = [Os t1.bandwidth(:,t1nodes) t2.bandwidth(:,t2nodes) ...
                  Os t1.bandwidth(:,t1leaves) t2.bandwidth(:,t2leaves) ...
                  Os t1.bandwidth(:,t1nodes+2*t1N*t1varBWs) ...
                  t2.bandwidth(:,t2nodes+2*t2.N*t2varBWs) Os ...
                  t1.bandwidth(:,t1leaves+2*t1N*t1varBWs) ...
                  t2.bandwidth(:,t2leaves+2*t2.N*t2varBWs) Os ...
                  t1.bandwidth(:,t1nodes+4*t1N*t1varBWs) ...
                  t2.bandwidth(:,t2nodes+4*t2.N*t2varBWs) Os ...
                  t1.bandwidth(:,t1leaves+4*t1N*t1varBWs) ...
                  t2.bandwidth(:,t2leaves+4*t2.N*t2varBWs) ];
else
  t1.bandwidth = [Os t1.bandwidth(:,t1nodes) t2.bandwidth(:,t2nodes) ...
                  Os t1.bandwidth(:,t1leaves) t2.bandwidth(:,t2leaves)];
end

%%%% Do stuff from calcStats because calcStats is protected.  Don't
% change calc stats or this won't work.

ax = max(t1.centers(:,t1root)+t1.ranges(:,t1root), ...
         t1.centers(:,t2root)+t1.ranges(:,t2root));
in = min(t1.centers(:,t1root)-t1.ranges(:,t1root), ...
         t1.centers(:,t2root)-t1.ranges(:,t2root));
t1.centers(:,1) = (ax+in)/2;
t1.ranges(:,1) = (ax-in)/2;

%calcuate weight
t1.weights(1) = t1.weights(t1root) + t1.weights(t2root);
W = sum(t1.weights(t1.N+1:2*t1.N));
t1.weights = t1.weights / W;  % normalize
t1w = t1.weights(t1root) / t1.weights(1);
t2w = t1.weights(t2root) / t1.weights(1);

%calculate mean
t1.means(:,1) = t1w*t1.means(:,t1root) + t2w*t1.means(:,t2root);

%calculate bandwidth
type = getType(t1);
if(strcmp(type, 'Gaussian'))
  t1.bandwidth(:,1) = t1w * (t1.bandwidth(:,t1root) + t1.means(:,t1root).^2) ...
      + t2w * (t1.bandwidth(:,t2root) + t1.means(:,t2root).^2) - ...
      t1.means(:,1).^2;
elseif(strcmp(type, 'Epanetchnikov'))
  t1.bandwidth(:,1) = sqrt(.5 * (t1w * (2*t1.bandwidth(:,t1root).^2 ...
                                                 + t1.means(:,t1root).^2) ...
                                + t2w * (2*t1.bandwidth(:,t2root).^2 ...
                                                   + t1.means(:,t2root).^2) ...
                                - t1.means(:,1).^2));
elseif(strcmp(type, 'Laplacian'))
  t1.bandwidth(:,1) = sqrt(5 * (t1w * (.2*t1.bandwidth(:,t1root).^2 ...
                                                + t1.means(:,t1root).^2) ...
                               + t2w * (.2*t1.bandwidth(:,t2root).^2 ...
                                                  + t1.means(:,t2root).^2) ...
                               - t1.means(:,1).^2));
else
  error(['unknown kernel type: ' type])
end

% take care of max and min BWs for variable bandwidths
if(varBWs)
  t1.bandwidth(:,1+2*t1.N) = max(t1.bandwidth(:,t1root+2*t1.N), t1.bandwidth(:,t2root+2*t1.N));
  t1.bandwidth(:,1+4*t1.N) = min(t1.bandwidth(:,t1root+4*t1.N), t1.bandwidth(:,t2root+4*t1.N));
end

t1n = 1;
t2n = 1 + length(t1nodes);
t1l = 1 + length(t2nodes);
t2l = 1 + length(t1nodes) + length(t1leaves);

% arrays are zero indexed
t1.lower = [t1.N addUints(t1.lower(t1nodes),t1l) ...
           addUints(t2.lower(t2nodes),t2l) 0 ...
           addUints(t1.lower(t1leaves),t1l) ...
           addUints(t2.lower(t2leaves),t2l)];
t1.upper = [2*t1.N-1 addUints(t1.upper(t1nodes),t1l) ...
           addUints(t2.upper(t2nodes),t2l) 0 ...
           addUints(t1.upper(t1leaves),t1l) ...
           addUints(t2.upper(t2leaves),t2l)];

t1leftch = addUints(t1.leftch, (t1.leftch < t1N) * t1n + ...
                     (t1.leftch >= t1N) * t1l);
t2leftch = addUints(t2.leftch, (t2.leftch < t2.N) * t2n + ...
                     (t2.leftch >= t2.N) * t2l);
t1rightch = addUints(t1.rightch, (t1.rightch < t1N) * t1n + ...
                     (t1.rightch >= t1N) .* (t1.rightch < 4e9) * t1l);
t2rightch = addUints(t2.rightch, (t2.rightch < t2.N) * t2n + ...
                     (t2.rightch >= t2.N) .* (t2.rightch < 4e9) * t2l);

t1.leftch = [d2uint(t1root) t1leftch(t1nodes) t2leftch(t2nodes) ...
             0 t1leftch(t1leaves) t2leftch(t2leaves)];
t1.rightch = [d2uint(t2root) t1rightch(t1nodes), t2rightch(t2nodes) ...
              0 t1rightch(t1leaves) t2rightch(t2leaves)];
t1.perm = [0 t1.perm(t1nodes) t2.perm(t2nodes) 0 t1.perm(t1leaves) ...
          addUints(t2.perm(t2leaves),length(t1leaves))];

function c = addUints(a, b)
c = uint32(double(a) + double(b));

function u = d2uint(d)
u = uint32(d - 1);

author: ShapeNet

repo: RenderForCNN


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

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

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

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

J = ((-1/m)*((y'*log(sigmoid(X*theta))) + ((1-y)'*log(1 - sigmoid(X*theta))))) + ((lambda/(2*m))*((theta' * theta)-(theta(1)^2)));

grad = ((1/m)*(X' * (sigmoid(X*theta) - y)))+((lambda/m)*(theta));

grad(1) = grad(1) - ((lambda/m)*(theta(1)));

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

end

author: UtkarshPathrabe


text:
function [h]=disp_console(h,m)
%
%
%

mm = get(h,'string');
mm{end+1} = m;
set(h,'string',mm);


text:
function neqn = p05_neqn ( )

%*****************************************************************************80
%
%% P05_NEQN returns the number of equations for problem P05.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NEQN, the number of equations.
%
  neqn = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% pdfbcalc_imagekld.m
% written by: Duncan Po
% Date: November 16, 2002
% calculate the Kublick-Liebler distance between two images based on
% on their hidden markov tree parameters
% Usage:    kld = pdfbcalc_imagekld(savename, sdir, model1, model2)
%           kld = pdfbcalc_imagekld(savename, sdir, image1, model2, imformat1)
%           kld = pdfbcalc_imagekld(savename, sdir, image1, image2, imformat1,
%                       imformat2)
% Inputs:   model1      - first Hidden Markov Tree Model
%           model2      - second Hidden Markov Tree Model
%           image1      - first image file name
%           image2      - second image file name
%           imformat1   - format of the first image file
%           imformat2   - format of the second image file
%           savename    - filename to save the model of image1
%                       - set to '' if no save desired
%           sdir        - full pathname to save the model 1 (also set to ''
%                           if savename == '')
% Output:   kld         - the Kullback Liebler distance between the trees
%           kld2        - KLD estimation using Monte Carlo method

%function [kld, kld2] = pdfbcalc_imagekld(savename, sdir, model1, model2, ...
%    imf1, imf2)
function [kld, kld2] = pdfbcalc_imagekld(savename, sdir, model1, model2, ...
    imf1, imf2)

% initialize the kld
kld = 0;
kld2 = 0;

if nargin == 6
    % this case is when the images are provided but not the models
    image1 = model1;
    image2 = model2;

    % train the HMT models first
    [model1, dummy] = pdfbtrainimagethmt(image1, imf1, '', 0.01);
    [model2, dummy] = pdfbtrainimagethmt(image2, imf2, '', 0.01);

elseif nargin == 5
    % this case is when we have one model and one image, need to
    % train a second model using the image provided
    image1 = model1;

    % train the HMT model
    [model1, dummy] = pdfbtrainimagethmt(image1, imf1, '', 0.01);

    %there is also the case for two models and no image when nargin==3
    %in that case, no training is needed
end;

if strcmp(savename, '') ~= 1
    file = sprintf('%s%s', sdir, savename);
    save(file,'model1');
end;

N = 60;
nlevel = model1{1}.nlevels;

for k = 1:nlevel
    levndir(k) = log2(length(model1{1}.stdv{k}));
end;

% treat each tree as independent and using chain rule of KLD for
% independent variables, simply add up the distances
for index = 1:length(model1)
    minkld = 10000;
    dbmodels = pdfbcreate_equiv_models(model1{index});
    for mmm = 1:length(dbmodels)
        tkld = pdfbcalc_KLD(nlevel, levndir, model2{index}, dbmodels{mmm});
        if tkld < minkld
            minkld = tkld;
        end;
    end;
    kld = kld + minkld;
end;

% Also apply the Monte Carlo method to verify the KLD estimation
for index = 1:length(model1)
    kld2 = kld2 + pdfbest_KLD(nlevel, levndir ,model2{index}, ...
        model1{index}, N);
end;

repo: mcbench-benchmarks


text:
function chrpak_test006 ( )

%*****************************************************************************80
%
%% TEST006 tests BINARY_TO_I4 and I4_TO_BINARY.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 June 2007
%
%  Author:
%
%    John Burkardt
%
  test_num = 4;

  i4_test = [ 21, -32, 2, 128 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST006\n' );
  fprintf ( 1, '  BINARY_TO_I4 converts a binary to an integer.\n' );
  fprintf ( 1, '  I4_TO_BINARY converts an integer to binary,\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I4 ==> BINARY ==> I4\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    i4 = i4_test(test);
    s = i4_to_binary ( i4 );
    j4 = binary_to_i4 ( s );

    fprintf ( 1, '  %8d  %8s  %8d\n', i4, s, j4 );

  end

  return
end


text:
function [FisoInd]=faceIsotropyIndex(F,V)
%% function for computing  triangular faces isotropy index.
% This is an index of the  "regularity" of the triangle. Equilateral triangle (perfectly regular)
%
% INPUTS:
% * F: nFaces-by-3 array representing the list of vertices for triangular faces
% * V: nVertices-by-3 array representing the 3d positions of the vertices of F in the reference positions
%
% OUTPUTS:
% * FisoInd: isotropy index for each face. value of 1 means perfectly
% equilateral. value of 0 means vertices are aligned (on a line)
%
%From the paper: Surface-Marker Cluster Design Criteria for 3-D Bone Movement Reconstruction (1997)
% Aurelio Cappozzo, Angelo Cappello, Ugo Della Croce, and Francesco Pensalfini
%%

FisoInd=zeros(size(F,1),1);
for iface=1:size(F,1)

    x1=V(F(iface,1),:)';
    x2=V(F(iface,2),:)';
    x3=V(F(iface,3),:)';

    if any(any(isnan([x1 x2 x3]))) % if any nan, iso index=nan
        FisoInd(iface)=NaN;
    else
        xa=(x1+x2+x3)/3;
        X=[x1-xa,x2-xa,x3-xa]; %cluster position model
        K=X*X'/3;
        Keig=real(eig(K)); %eigenvalues of X
        FisoInd(iface)=2*Keig(2)/(Keig(2)+Keig(3)); % isotropy index
    end
end

end

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

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC