text:
function bdt_analysis() %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Workspace initialization
%  The workspace is cleaned and some parameters and deffined
clear all ; % this removes all variables stored in your current workspace

% PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The next parameters are choosen to produce an optimal solution. I would
% recommend not to modify them, but they could be changed by the user to
% produce different scnearios.
%
% (1) <NoDT> The Number of Decission Trees to construct the embedded model.
%     For example:
%     NoDT = 100 ;
      NoDT = 100 ;
%
% (2) <seed> This parameter is the seed of the random-number generator. By
%     fixing it, we have controlled the randomness of the results and we
%     could reproduce it at any moment. The seed could have any integer
%     value, but I would recommend to reset it always by ussing zero.
%     For example:
%     seed = 0 ;
      seed = 0 ;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all force ; % this closes all plotted figures
clc ; % this clears the command window log
rng(seed) ; % random seed is set here

disp( ['==============================================='] ) ;
disp( ['= Classification using Bagged Decission Trees ='] ) ;
disp( ['==============================================='] ) ; disp( char(10) ) ;

%% Loading and examinating the data
%  Data are read from an Excel or CSV file (genes at columns and samples at rows)
%  which must have the next structure:
%  - First column: the name of the class (cell-type, treatment, stage, etc)
%  - Second column: the name of the sample
%  - First row: a header with the name of the genes at each column. Notice
%  that the first two cells of the header corresponds to the sample's
%  class and samples's name.
tic ;  disp( [' - Loading data...'] ) ; disp( char(10) ) ;
[ input_file,input_path ] = uigetfile( {'*.csv','Comma Sparated Values (*.csv)';'*.xls','Excel 97-2003 (*.xls)';'*.xlsx','Excel 2010 (*.xlsx)'},'MultiSelect','off' ) ;
[ vector ] = strsplit( input_file,'.' ) ;
file_name = vector(1) ;
extension = vector(2) ;
if ( strcmp(extension,'csv') )
  T = readtable( [input_path,input_file],'Delimiter',',','ReadVariableNames',false ) ;
  classes_names = T{2:end,1} ;
  samples = T{2:end,2} ;
  if ( sum(sum(isnan(char(T{:,end})))) > 0 )
    genes = T{1,3:end-1} ;
    data = str2double(T{2:end,3:end-1}) ;
  else
    genes = T{1,3:end} ;
    data = str2double(T{2:end,3:end}) ;
  end%if
  clear T ;
else
  if ( strcmp(extension,'xls') ||  strcmp(extension,'xlsx') )
    [ data,txt ] = xlsread( [input_path,input_file] ) ;
    classes_names = txt(2:end,1) ;
    samples = txt(2:end,2) ;
    genes = txt(1,3:end) ;
    clear txt ;
  else
    disp('   file-format error.')
    return ;
  end%if
end%if
whos('data','samples','classes_names','genes') ;
%  "data" is the matrix containing the normalized expression profiles from
%  single-cell qPCR data, with genes as columns and samples as rows. The
%  name of the genes at each column are in "genes", a row-vector.
%  Similarly, the name of each sample at each row of the "data" matrix are
%  stored in "cell_names", a column-vector. The names on these vectors must
%  be unique. The name of the classes at each row are stored in
%  "classes_names". This column vector has non-unique classes names.
N = size(data,1) ;
G = size(data,2) ;
classes_unique = unique( sort(classes_names) ) ;
K = size(classes_unique(:),1) ;
classes_keys = zeros(N,1) ;
classes_matrix = zeros(N,N) ;
% The population of classes can be introduced as a prior knowledge by
% indicating the frecuency of each class. A non informative prior should be
% a constant value for all the classes, for instances:
  classes_prior = ones(K,1) ;
% Alternatively you can use your custom values.
% classes_prior = [ 1.42 , 1.42 , 0.4 , 0.25 , 0.5 , 0.20 , 0.07 , 1.42 , 1.42 ]' ;
  for k = 1:K
    indexes = strcmp( classes_unique(k),classes_names ) ;
    classes_keys(indexes) = k ;
    classes_matrix = classes_matrix + 1.0*indexes*indexes' ;
  end%for
%  "classes_keys" is a column vector in which the classes names are mapped
%  to a key column-vector of integers.
fprintf('   dataset with %i samples, %i genes and %i unique classes.',N,G,K ) ; disp( char(10) ) ;

%% Building the model
%  The Bagged Decission Tree builds an ensemble with a fixed Number Of
%  Deicission Trees (NODT). The ensemble is generated for the whole dataset and the Out-Of-Bag (OOB)
%  infromation is stored for evaluation purposes.
tic ;  disp( [' - Building the model with ',num2str(NoDT),' trees...'] ) ;
dbt = TreeBagger( NoDT,data,classes_names,'names',genes,'method','classification','NVarToSample','all','oobvarimp','on','oobpred','on','Prior',classes_prior ) ;
fprintf('   ( +%0.2f s ) done!',toc) ; disp( char(10) ) ;

%% Biased evaluation of the model
%  This is a first test of the model, consisting of the prediction the
%  classes of the same dataset used for the trainning. This is a biased
%  evaluation that is not significant.
[ predicted_classes predicted_scores ] = predict( dbt,data ) ;
%  First, we plot the scores of all the predicted classes for each sample.
%  The color key corresponds to the actual class of each sample.
figure(1) ;  set(1,'WindowStyle','docked') ; clf ;
bp = nan(size(predicted_scores)) ;
colors = hsv(K+1) ;
h = [] ;
for k = 1:K
  indexes = classes_keys == k ;
  bp(indexes,k) = predicted_scores(indexes,k) ;
end%for
subplot(2,1,1) ; cla ;
hold on ;
  boxplot(bp,'Colors',colors) ;
hold off ;
grid on ;
xlim([1 K]) ; xlabel( 'Classes' ) ;
ylim([0 1]) ; ylabel( 'Scores' ) ;
title('Scores of predicted classes (biased evaluation).') ;
set(gca,'XTick',[1:K]) ;
set(gca,'XTickLabel',classes_unique) ;
%  Second, we compute the confusion matrix (CM). The main diagonal of the
%  CM holds the number of correct classified samples. Otherwise, the
%  elements out of the main diagonal denotes wrongly classified samples.
[ cm cm_labels ] = confusionmat( classes_names,predicted_classes,'order',classes_unique ) ;
%  For a full inspection of the matrix, use the next command.
disp( dataset({cm,cm_labels{:}},'obsnames',cm_labels) ) ;
%  Alternatively, use a stacked-bar plot to visualize the proportions of
%  predicted classes for each class.
subplot(2,1,2) ; cla ;
hold on ;
  bar( cm,'stacked' ) ;
hold off ;
grid on ;
colormap( colors ) ;
xlabel( 'Classes' ) ; xlim([0 K]+0.5) ;
ylabel( 'Number of samples' ) ;
title('Fraction of predicted classes versus the actual ones (biased evaluation).') ;
legend( cm_labels,'EdgeColor',[1 1 1],'Location','NO','Orientation','horizontal' ) ;
set(gca,'XTick',[1:K]) ;
set(gca,'XTickLabel',cm_labels) ;

%% Evaluating the model with the OOB samples
%  This evaluation takes into account the out-of-bag (OOB) samples, that
%  have been not considered to build the model.
[ predicted_classes predicted_scores ] = oobPredict( dbt ) ;
figure(2) ; set(2,'WindowStyle','docked') ; cla ;
bp = nan(size(predicted_scores)) ;
colors = hsv(K+1) ;
h = [] ;
hold on ;
  for k = 1:K
    indexes = classes_keys == k ;
    bp(indexes,k) = predicted_scores(indexes,k) ;
  end%for
subplot(2,1,1) ;
hold on ;
  boxplot(bp,'Colors',colors) ;
hold off ;
grid on ;
xlim([1 K]) ; xlabel( 'Classes' ) ;
ylim([0 1]) ; ylabel( 'Scores' ) ;
title('Scores of predicted classes (OOB evaluation).') ;
set(gca,'XTick',[1:K]) ;
set(gca,'XTickLabel',classes_unique) ;
[ cm cm_labels ] = confusionmat( classes_names,predicted_classes,'order',classes_unique ) ;
disp( dataset({cm,cm_labels{:}},'obsnames',cm_labels) ) ;
subplot(2,1,2) ; cla ;
hold on ;
  bar( cm,'stacked' ) ;
hold off ;
grid on ;
colormap( colors ) ;
xlabel( 'Classes' ) ; xlim([0 K]+0.5) ;
ylabel( 'Number of samples' ) ;
title('Fraction of predicted classes versus the actual ones (OOB evaluation).') ;
legend( cm_labels,'EdgeColor',[1 1 1],'Location','NO','Orientation','horizontal' ) ;
set(gca,'XTick',[1:K]) ;
set(gca,'XTickLabel',cm_labels) ;

%% Examining the errors of the individual models
%  To test the model objectively we use the OOB error (its missclasication
%  probability) and its margin (the difference between the score for the
%  actual class and the largest one for other classes) for each
%  individual tree and theirs cumulative values for the emsembled model.
figure(3) ; set(3,'WindowStyle','docked') ; cla ;
subplot(2,1,1) ; cla ;
hold on ;
  plot( oobError(dbt,'mode','individual'),'ro','Color',[1 0.2 0] ) ;
  curve=smooth(oobError(dbt,'mode','cumulative')) ;
  plot( curve,'r-','Color',[1 0.2 0] ) ;
  index=find(abs(diff(curve))<0.001*max(curve(:)),1) ;
  plot( [1,NoDT],curve(index)*[1 1],'r--','Color',[1 0.2 0] ) ;
hold off ;
grid on ;
xlabel( 'Number of trees in the ensemble' ) ; xlim([1 NoDT]) ;
ylabel( 'Misclassification probability' ) ; ylim([0 1]) ;
title( 'Missclassification probability for the OOB samples using an incremental ensembling.' ) ;
legend( {'Individual';'Cumulative';'Steady level'},'EdgeColor',[1 1 1],'Location','NE','Orientation','horizontal' ) ;
%  The margin can be considered as a confidence interval of the score. The
%  highest the margin, the bigger the confidence on the predicted class.
subplot(2,1,2) ; cla ;
hold on ;
  plot( oobMeanMargin(dbt,'mode','individual'),'bo','Color',[0 0.6 1] ) ;
  curve=smooth(oobMeanMargin(dbt,'mode','cumulative')) ;
  plot( curve,'r-','Color',[0 0.6 1] ) ;
  index=find(abs(diff(curve))<0.001*max(curve(:)),1) ;
  plot( [1,NoDT],curve(index)*[1 1],'r--','Color',[0 0.6 1] ) ;
hold off ;
grid on ;
xlabel( 'Number of trees in the ensemble' ) ; xlim([1 NoDT]) ;
ylabel( 'Margin' ) ; ylim([0 1]) ;
title( 'Mean of the margin for the OOB samples using an incremental ensembling.' ) ;
legend( {'Individual';'Cumulative';'Steady level'},'EdgeColor',[1 1 1],'Location','NE','Orientation','horizontal' ) ;

%% Examining the correlation among the predicted classes
%  The proximitry matrix is defined as the fraction of trees in the ensemble
%  for which any two observations land on the same leaf.
dbt = fillProximities( dbt ) ;
hold on ;
  HeatMap( dbt.Proximity,'ColumnLabels',classes_names,'RowLabels',classes_names,'Colormap',redgreencmap(20) ) ;
hold off ;
title('Proximity matrix.') ;

%% Examining the confidence of the invidual variables
%  The relevance of each gene on the emsenble may be computed as the
%  predictor importance. The predictor importance averages the changes in
%  the risk due to split on every predictor at each node for the whole
%  ensemble. This risk depends on the seleted split criterion (typically
%  GDI that measure the impurity of the node). "DeltaCritDecisionSplit"
%  measure the changes in the split criterion summed over splits on each
%  variable, averaged across the entire ensemble.
[ sorted_improvement_criteria sorted_improvement_keys ] = sort( dbt.DeltaCritDecisionSplit,'descend' ) ;
figure(4) ; set(4,'WindowStyle','docked') ; cla ;
hold on ;
  bar( sorted_improvement_criteria,'FaceColor',[0 0.6 1],'EdgeColor',[0 0.4 0.6] ) ;
hold off ;
grid on ;
xlabel( 'Genes' ) ; xlim([0 G]+0.5) ;
ylabel( 'Changes in the split criterion' ) ;
title('Predictor improvement in the split criterion (descending sorted).') ;
set( gca,'XTick',[1:G] ) ;
set( gca,'XTickLabel',genes(sorted_improvement_keys) ) ;

%% Examining the errors of the individual variables
%  Alternaitvely, the relevance of each variable may be computed by
%  "OOBPermutedVarDeltaError". This variable measures the increase in the
%  OOB prediction error if the values of that variables are permutted for
%  every tree, then averaged over the entire ensemble and lastly divided by
%  the standard deviation.
[ sorted_error_criteria sorted_error_keys ] = sort( dbt.OOBPermutedVarDeltaError,'descend' ) ;
figure(5) ; set(5,'WindowStyle','docked') ; cla ;
hold on ;
  bar( sorted_error_criteria,'FaceColor',[0 0.6 1],'EdgeColor',[0 0.4 0.8] ) ;
hold off ;
grid on ;
xlabel( 'Genes' ) ; xlim([0 G]+0.5) ;
ylabel( 'Effect on the prediction error' ) ;
title( 'Predictor improvement in the OOB error (descending sorted).' ) ;
legend( {'Error rate'},'EdgeColor',[1 1 1],'Location','NE' ) ;
set(gca,'XTick',[1:G]) ;
set(gca,'XTickLabel',genes(sorted_error_keys)) ;

%% Plotting the performance curves
figure(6) ; set(6,'WindowStyle','docked') ; cla ;
text_legend = [] ;
hold on ;
for k = 1:K
  [ fpr tpr threshold AUC(k) ] = perfcurve( classes_names,predicted_scores(:,k),classes_unique{k} ) ;
  plot(fpr,tpr,'Color',colors(k,:)) ;
  text_legend = [ text_legend ; {[classes_unique{k},' (',num2str(AUC(k)),')']} ] ;
end%for
 plot([0 1],[0 1],'g--','Color',[1 1 1]*0.5)
hold off ;
grid on ;
legend( [ text_legend(:) ; 'ND line' ] ,'EdgeColor',[1 1 1],'Location','SE' ) ;
title('ROC curve (TPR and FPR trade off).') ;
xlabel( '1-Specificity (FPR)' ) ; xlim([0 1]) ;
ylabel( 'Sensitivity (TPR)' ) ; ylim([0 1]) ;
figure(7) ; set(7,'WindowStyle','docked') ; cla ;
text_legend = [] ;
hold on ;
for k = 1:K
  [ ppv tpr threshold AUC(k) ] = perfcurve( classes_names,predicted_scores(:,k),classes_unique{k},'xcrit','ppv' ) ;
  plot(tpr,ppv,'Color',colors(k,:)) ;
  text_legend = [ text_legend ; {[classes_unique{k},' (',num2str(AUC(k)),')']} ] ;
end%for
hold off ;
grid on ;
legend( [ text_legend(:) ] ,'EdgeColor',[1 1 1],'Location','SW' ) ;
title('PR curve (PPV and TPR trade off).') ;
xlabel( 'Recall (TPR)' ) ; xlim([0 1]) ;
ylabel( 'Precission (PPV)' ) ; ylim([0 1]) ;

%% Plotting the tree and loading super cells
% Plotting the embedded tree
view(dbt.Trees{:},'mode','graph')
% Tracking a subset of cells along the tree
choice = questdlg(['Do you want to track a superclass along the tree?'],'Track Superclass','Yes','No','Yes');
switch choice
case 'Yes'
%  Data should be provided in a column TXT file with no header.
tic ;  disp( [' - Loading data...'] ) ; disp( char(10) ) ;
[ input_file_super,input_path_super ] = uigetfile( {'*.txt','One column text format (*.txt)'},'MultiSelect','off','Select File to Open',input_path ) ;
super_samples = readtable( [input_path_super,input_file_super],'Delimiter',',','ReadVariableNames',false ) ;
super_samples = super_samples{:,1} ;
indexes = zeros(N,1) ;
for k = 1:size(super_samples(:),1)
  indexes = indexes + strcmp( samples,super_samples(k) ) ;
end%for
super_indexes = indexes == 1 ;
import bioma.data.*
dm = DataMatrix( data(super_indexes,:),samples(super_indexes),genes ) ;
super_classes = classes_names(super_indexes) ;
super_classes_uinque = unique(sort(super_classes)) ;
colors = hsv( size(super_classes_uinque(:),1) ) ;
markers = { 'o' ; 's' ; 't' ; 'x' ; 'h' ; 'p' }' ;
%dm = DataMatrix(data,classes_names,genes) ;
embedded = dbt.Trees{:} ;
node = 1 ;
leave = 1 ;
splitting_variables = {} ;
splitting_values = [] ;
parents = ones(size(dm,1),1) ;
tree = [ 0 ] ;
family = parents(:)' ;
[ tree , family , splitting_variables , splitting_values ] = growtree( dm,embedded,node,leave,parents,tree,family,splitting_variables,splitting_values ) ;
tree = tree(:)' ;
count = size(tree,2) ;
[ x , y ] = treelayout( tree ) ;
figure(8) ; cla ; hold on ; xlim([0 1]) ; ylim([0 1]) ;
    h=[] ; legend_text = [] ;
    for k = 1:size(super_classes_uinque(:),1)
      legend_text = [ legend_text ; super_classes_uinque(k) ] ;
      h(k) = plot( 2,2,'k.-','Color',colors(k,:) ) ;
    end%for
for k = 1:1:size(splitting_variables(:),1)
    if sum( tree == k ) > 0
      indexes2 = find( tree==k ) ;
      for j=1:size(indexes2(:),1)
        plot([x(k) x(indexes2(j))],[y(k) y(indexes2(j))],'k--','Color',[1 1 1]*0.5) ;
      end%if
      if sum( tree == k ) == 2
        text( -0.01+x(k),0.02+y(k),splitting_variables(k) ) ;
        population = [] ;
        for i=1:size(super_classes_uinque(:),1)
          population = [ population ; sum( family(k,:).*strcmp(super_classes,super_classes_uinque(i))' ) ] ;
        end%for
        for i=1:size(population(:),1)
          if sum( population(i) ) > 0
            plot( x(k)+0.005*i-0.005*size(super_classes_uinque(:),1)/1.75,y(k),'ro--','Color',colors(i,:),'MarkerSize',50*population(i)/size(super_classes,1) ) ;
          end%if
        end%for
      end%if
    else
      %plot(x(k),y(k),'r.') ;
        population = [] ;
        for i=1:size(super_classes_uinque(:),1)
          population = [ population ; sum( family(k,:).*strcmp(super_classes,super_classes_uinque(i))' ) ] ;
          if population(i) > 0
            plot( x(k),y(k)+0.005*i-0.005*size(super_classes_uinque(:),1)/1.75,'ro--','Color',colors(i,:),'MarkerSize',50*population(i)/size(super_classes,1) ) ;
          end%if
        end%for
    end%if
    drawnow ;
end%for
hold off ;
legend( h(:),legend_text,'EdgeColor',[1 1 1],'Location','NE' ) ; set(gca,'XTick',[]) ; set(gca,'YTick',[]) ;
end%switch

set(8,'WindowStyle','docked') ;

% http://stackoverflow.com/questions/5065051/add-node-numbers-get-node-locations-from-matlabs-treeplot

end%function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ tree , family , splitting_variables , splitting_values ] = growtree( dm,embedded,node,leave,parents,tree,family,splitting_variables,splitting_values )

  if ~( strcmp(embedded.CutVar(leave),'') )
    c1 = ( parents.*( dm( :,embedded.CutVar(leave) ) < embedded.CutPoint(leave) ) ) == 1 ;
    node1 = 0 ;
    c2 = ( parents.*( c1 == 0 ) ) == 1 ;
    node2 = 0 ;
    if sum( c1 ) >= 2 % minimun number of children at last leave
      tree = [ tree(:) ; node ] ;
      family = [ family ; c1(:)' ] ;
      node1 = size(tree(:),1) ;
      splitting_variables = [ splitting_variables ; embedded.CutVar(leave) ] ;
      splitting_values = [ splitting_values ; embedded.CutPoint(leave) ] ;
      [ tree , family , splitting_variables , splitting_values ] = growtree( dm,embedded,node1,embedded.Children(leave,1),c1,tree,family,splitting_variables,splitting_values ) ;
    end%if
    if sum( c2 ) >= 2 % minimun number of children at last leave
      tree = [ tree(:) ; node ] ;
      family = [ family ; c2(:)' ] ;
      node2 = size(tree(:),1) ;
      splitting_variables = [ splitting_variables ; embedded.CutVar(leave) ] ;
      splitting_values = [ splitting_values ; embedded.CutPoint(leave) ] ;
      [ tree , family , splitting_variables , splitting_values ] = growtree( dm,embedded,node2,embedded.Children(leave,2),c2,tree,family,splitting_variables,splitting_values ) ;
    end%if
  end%if

end%function

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function PLOT_Stop_Action(States,Params,Trajectory)

% This function is used to plot the "stop action" version of the simulation.
% This is analgous to running the simulation slowly with the "hold on"
% feature set.

Fig_Num = 4;
N_Frames = 20;

Npts = Params.Traj.Npts;

x = States(1,:);
y = States(2,:);
th = States(3,:);
phi = States(4,:);

Lt = Params.Dyn.Lt;
Lc = Params.Dyn.Lc;

Nx = length(x);

Display = linspace(1,Nx,N_Frames);
Display = round(Display);

%[dStates, A, B, Positions] = Derive_EoM()

Pa = [...   %Position of the back of the trailer
    x;
    y];

Pb = [...    %Position of the back of the cab
 x - Lt*sin(th);
 y + Lt*cos(th)];

Pc = [...   %Position of the front of the cab
 x - Lc*(cos(phi).*sin(th) + cos(th).*sin(phi)) - Lt*sin(th);
 y + Lc*(cos(phi).*cos(th) - sin(phi).*sin(th)) + Lt*cos(th)];

if strcmp(Params.Sim.Direction,'reverse')
    Dir = 2;
else
    Dir = 1;
end
Road = zeros(2,Npts);
for i=1:Npts
    Road(:,i) = Trajectory.States{Dir,i}(1:2);
end

All_Points = [Pa,Pb,Pc,Road];

Xmin = min(All_Points(1,:));
Xmax = max(All_Points(1,:));
Ymin = min(All_Points(2,:));
Ymax = max(All_Points(2,:));

figure(Fig_Num)
clf;
hold on
axis([Xmin,Xmax,Ymin,Ymax])
axis equal
title('Simulation of the Tractor Trailer Truck')
plot(Road(1,:),Road(2,:),'k','LineWidth',Params.Sim.RoadWidth)

for i=Display
    Trailer_X = [Pa(1,i); Pb(1,i)];
    Trailer_Y = [Pa(2,i); Pb(2,i)];
    Cab_X = [Pb(1,i); Pc(1,i)];
    Cab_Y = [Pb(2,i); Pc(2,i)];

    plot(Trailer_X,Trailer_Y,'b','LineWidth',20)
    plot(Cab_X, Cab_Y,'g','LineWidth',17)
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function [data] = ft_resampledata(cfg, data)

% FT_RESAMPLEDATA performs a resampling or downsampling of the data
%
% Use as
%   [data] = ft_resampledata(cfg, data)
%
% The data should be organised in a structure as obtained from the FT_PREPROCESSING
% function. The configuration should contain
%   cfg.resamplefs      = frequency at which the data will be resampled
%   cfg.detrend         = 'no' or 'yes', detrend the data prior to resampling (no default specified, see below)
%   cfg.demean          = 'no' or 'yes', whether to apply baseline correction (default = 'no')
%   cfg.baselinewindow  = [begin end] in seconds, the default is the complete trial (default = 'all')
%   cfg.feedback        = 'no', 'text', 'textbar', 'gui' (default = 'text')
%   cfg.trials          = 'all' or a selection given as a 1xN vector (default = 'all')
%   cfg.sampleindex     = 'no' or 'yes', add a channel with the original sample indices (default = 'no')
%
% Instead of specifying cfg.resamplefs, you can also specify a time axis on which you
% want the data to be resampled. This is useful for merging data from two acquisition
% devices, after resampledata you can call FT_APPENDDATA to concatenate the channels
% from the different acquisition devices.
%   cfg.time        = cell-array with one time axis per trial (i.e. from another dataset)
%   cfg.method      = interpolation method, see INTERP1 (default = 'pchip')
%   cfg.extrapval   = extrapolation behaviour, scalar value or 'extrap' (default = as in INTERP1)
%
% When you specify cfg.method='resample' an implicit anti-aliasing low pass filter is applied prior
% to the resampling. You can also explicitly specify an anti-aliasing low pass filter. This is adviced
% when downsampling using any other method than 'resample', but also when strong noise components are
% present just above the new Nyquist frequency.
%   cfg.lpfilter    = 'yes' or 'no' (default = 'no')
%   cfg.lpfreq      = scalar value for low pass frequency (there is no default, so needs to be always specified)
%   cfg.lpfilttype  = string, filter type (default is set in ft_preproc_lowpassfilter)
%   cfg.lpfiltord   = scalar, filter order (default is set in ft_preproc_lowpassfilter)
%
% More documentation about anti-alias filtering can be found in this FAQ on the FieldTrip website.
%
% Previously this function used to detrend the data by default to avoid edge artifacts fue to the
% anti-aliassing filter. Detrending is fine for removing slow drifts in data prior to frequency analysis,
% but not recommended if you want to look at ERPs or ERFs. Therefore the old default value 'yes' has been
% removed; you now explicitly have to specify whether you want to detrend or not.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_PREPROCESSING, FT_APPENDDATA, FT_PREPROC_LOWPASSFILTER, ESAMPLE, DOWNSAMPLE, INTERP1

% Copyright (C) 2003-2006, FC Donders Centre, Markus Siegel
% Copyright (C) 2004-2022, FC Donders Centre, Robert Oostenveld
% Copyright (C) 2022, DCCN, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar data
ft_preamble provenance data

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% ft_checkdata is done further down

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden', {'trial'}); % prevent accidental typos, see issue 1729
cfg = ft_checkconfig(cfg, 'renamed', {'blc', 'demean'});
cfg = ft_checkconfig(cfg, 'renamed', {'resamplemethod', 'method'});
cfg = ft_checkconfig(cfg, 'renamed', {'fsample', 'resamplefs'});

% set the defaults
cfg.resamplefs       = ft_getopt(cfg, 'resamplefs',      []);
cfg.time             = ft_getopt(cfg, 'time',            {});
cfg.factor           = ft_getopt(cfg, 'factor',          {});
cfg.detrend          = ft_getopt(cfg, 'detrend',         'no');
cfg.demean           = ft_getopt(cfg, 'demean',          'no');
cfg.baselinewindow   = ft_getopt(cfg, 'baselinewindow',  'all');
cfg.feedback         = ft_getopt(cfg, 'feedback',        'text');
cfg.trials           = ft_getopt(cfg, 'trials',          'all', 1);
cfg.method           = ft_getopt(cfg, 'method',          []);
cfg.sampleindex      = ft_getopt(cfg, 'sampleindex',     'no');
cfg.extrapval        = ft_getopt(cfg, 'extrapval',       []);
cfg.lpfilter         = ft_getopt(cfg, 'lpfilter');

% store original datatype
convert = ft_datatype(data);

% check if the input data is valid for this function, this will convert it to raw if needed
data = ft_checkdata(data, 'datatype', {'raw+comp', 'raw'}, 'feedback', 'yes');

if isempty(cfg.method) && ~isempty(cfg.time)
  % see INTERP1, shape-preserving piecewise cubic interpolation
  cfg.method = 'pchip';
elseif isempty(cfg.method)
  % see RESAMPLE
  cfg.method = 'resample';
end

usefsample = any(strcmp(cfg.method, {'resample', 'downsample', 'decimate', 'mean', 'median'}));
usetime    = ~usefsample;

% select trials of interest
tmpcfg = keepfields(cfg, {'trials', 'showcallinfo', 'trackcallinfo', 'trackusage', 'trackdatainfo', 'trackmeminfo', 'tracktimeinfo', 'checksize'});
data   = ft_selectdata(tmpcfg, data);
% restore the provenance information
[cfg, data] = rollback_provenance(cfg, data);

if strcmp(cfg.sampleindex, 'yes') && isfield(data, 'sampleinfo')
  data.label{end+1} = 'sampleindex';
  for i=1:size(data.sampleinfo,1)
    % this works for one or more trials
    data.trial{i}(end+1,:) = data.sampleinfo(i,1):data.sampleinfo(i,2);
  end
elseif strcmp(cfg.sampleindex, 'yes')
  ft_warning('no sampleinfo present, cannot add sampleindex as channel');
end

% sampleinfo, if present, becomes invalid because of the resampling
if isfield(data, 'sampleinfo')
  data = rmfield(data, 'sampleinfo');
end

if usefsample && usetime
  ft_error('you should either specify cfg.resamplefs or cfg.time')
end

% remember the original sampling frequency in the configuration
cfg.origfs = double(data.fsample);

% set this to nan, it will be updated later on
data.fsample = nan;

if isempty(cfg.lpfilter), cfg.lpfilter = 'no'; end
dolpfilt = istrue(cfg.lpfilter);
if dolpfilt
  cfg.lpfilttype = ft_getopt(cfg, 'lpfilttype');
  cfg.lpfiltord  = ft_getopt(cfg, 'lpfiltord');
  cfg            = ft_checkconfig(cfg, 'required', 'lpfreq');
end

if usefsample
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % resample/downsample based on new sampling frequency
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  ntr = length(data.trial);
  nchan  = numel(data.label);

  ft_progress('init', cfg.feedback, 'resampling data');
  [fsorig, fsres] = rat(cfg.origfs./cfg.resamplefs); %account for non-integer fs
  cfg.resamplefs  = cfg.origfs.*(fsres./fsorig); %get new fs exact

  % make sure that the resampled time axes are aligned (this is to avoid rounding
  % errors in the time axes). this procedure relies on the fact that resample assumes
  % all data outside the data window to be zero anyway. therefore, padding with zeros
  % (to the left and right) before resampling does not hurt
  begsample = zeros(ntr, 1);
  endsample = zeros(ntr, 1);
  for itr = 1:ntr
    begsample(itr) = round(cfg.origfs * data.time{itr}(1));
    endsample(itr) = round(cfg.origfs * data.time{itr}(end));
  end
  begpad = begsample-min(begsample);
  endpad = max(endsample)-endsample;

  if any(begpad~=0) || any(endpad~=0)
    ft_warning('not all trials have the same time axis; data will be zero-padded prior to resampling to avoid rounding issues in the resampled time axes');
  end

  if any(strcmp(cfg.method, {'downsample', 'mean', 'median'}))
    ft_warning('using cfg.method = ''%s''; only use this if you have applied an anti-aliasing filter prior to downsampling!', cfg.method);
  end

  if any(strcmp(cfg.method, {'decimate', 'downsample', 'mean', 'median'}))
    if mod(fsorig, fsres) ~= 0
      ft_error('the new sampling rate needs to be an integer division of the original sampling rate');
    end
  end

  for itr = 1:ntr
    ft_progress(itr/ntr, 'resampling data in trial %d from %d\n', itr, ntr);

    olddat = data.trial{itr};
    oldtim = data.time{itr};

    % detrending is in general not recommended
    if istrue(cfg.detrend)
      if ~strcmp(cfg.baselinewindow, 'all')
        olddat = ft_preproc_detrend(olddat, nearest(oldtim, cfg.baselinewindow(1)), nearest(oldtim, cfg.baselinewindow(2)));
      else
        olddat = ft_preproc_detrend(olddat);
      end
    end

    % remove the mean to avoid edge effects when there's a strong offset, the cfg.demean option is dealt with below
    if ~strcmp(cfg.baselinewindow, 'all')
      [olddat, bsl] = ft_preproc_baselinecorrect(olddat, nearest(oldtim, cfg.baselinewindow(1)), nearest(oldtim, cfg.baselinewindow(2)));
    else
      [olddat, bsl] = ft_preproc_baselinecorrect(olddat);
    end

    if istrue(cfg.lpfilter)
      olddat = ft_preproc_lowpassfilter(olddat, cfg.origfs, cfg.lpfreq, cfg.lpfiltord, cfg.lpfilttype);
    end

    % pad the data with zeros on both sides
    olddat = [zeros(nchan, begpad(itr)) olddat zeros(nchan, endpad(itr))];
    oldtim = ((begsample(itr)-begpad(itr)):(endsample(itr)+endpad(itr))) / cfg.origfs;

    % perform the resampling
    if strcmp(cfg.method, 'downsample')
      if isa(olddat, 'single')
        % temporary convert this trial to double precision
        newdat = transpose(single(downsample(double(transpose(olddat)),fsorig/fsres)));
      else
        newdat = transpose(downsample(transpose(olddat),fsorig/fsres));
      end

    elseif strcmp(cfg.method, 'resample')
      if isa(olddat, 'single')
        % temporary convert this trial to double precision
        newdat = transpose(single(resample(double(transpose(olddat)),fsres,fsorig)));
      else
        newdat = transpose(resample(transpose(olddat),fsres,fsorig));
      end

    elseif strcmp(cfg.method, 'decimate')
      if isa(olddat, 'single')
        % temporary convert this trial to double precision
        newdat = transpose(single(my_decimate(double(transpose(olddat)),fsorig/fsres)));
      else
        newdat = transpose(my_decimate(transpose(olddat),fsorig/fsres));
      end

    elseif strcmp(cfg.method, 'mean')
      if isa(olddat, 'single')
        % temporary convert this trial to double precision
        newdat = transpose(single(my_mean(double(transpose(olddat)),fsorig/fsres)));
      else
        newdat = transpose(my_mean(transpose(olddat),fsorig/fsres));
      end

    elseif strcmp(cfg.method, 'median')
      if isa(olddat, 'single')
        % temporary convert this trial to double precision
        newdat = transpose(single(my_median(double(transpose(olddat)), fsorig/fsres)));
      else
        newdat = transpose(my_median(transpose(olddat), fsorig/fsres));
      end

    else
      ft_error('unknown method ''%s''', cfg.method);
    end

    % add back the mean
    if ~strcmp(cfg.demean, 'yes')
      nsmp   = size(newdat,2);
      newdat = newdat + bsl(:,ones(1,nsmp));
    end

    % compute the new time axis, assuming that it starts at the same time
    nsmp   = size(newdat,2);
    newtim = (0:(nsmp-1))/cfg.resamplefs;

    % the middle of the time bin represented by the first samples are not aligned
    % the new time axis can be shifted by a sub-sample amount
    shift  = mean(oldtim) - mean(newtim);
    newtim = newtim + shift;

    if begpad(itr)>0 || endpad(itr)>0
      % un-pad the data
      sel = (1+round(begpad(itr)*cfg.resamplefs/cfg.origfs)):(length(newtim)-round(endpad(itr)*cfg.resamplefs/cfg.origfs));
      newtim = newtim(   sel);
      newdat = newdat(:, sel);
    end

    data.time{itr}  = newtim;
    data.trial{itr} = newdat;

  end % for itr
  ft_progress('close');

  % specify the new sampling frequency in the output
  data.fsample = cfg.resamplefs;

elseif usetime
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % resample based on the specified new time axes for each trial
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if isempty(cfg.extrapval)
    if strcmp(cfg.method, 'spline') || strcmp(cfg.method, 'pchip')
      cfg.extrapval = 'extrap';
    else
      cfg.extrapval = nan;
    end
  end

  ntr = length(data.trial);

  ft_progress('init', cfg.feedback, 'resampling data');
  for itr = 1:ntr
    ft_progress(itr/ntr, 'resampling data in trial %d from %d\n', itr, ntr);

    olddat = data.trial{itr};
    oldtim = data.time{itr};

    % detrending is in general not recommended
    if istrue(cfg.detrend)
      if ~strcmp(cfg.baselinewindow, 'all')
        olddat = ft_preproc_detrend(olddat, nearest(oldtim, cfg.baselinewindow(1)), nearest(oldtim, cfg.baselinewindow(2)));
      else
        olddat = ft_preproc_detrend(olddat);
      end
    end

    % always remove the mean to avoid edge effects when there's a strong offset, the cfg.demean option is dealt with below
    if ~strcmp(cfg.baselinewindow, 'all')
      [olddat, bsl] = ft_preproc_baselinecorrect(olddat, nearest(oldtim, cfg.baselinewindow(1)), nearest(oldtim, cfg.baselinewindow(2)));
    else
      [olddat, bsl] = ft_preproc_baselinecorrect(olddat);
    end

    if istrue(cfg.lpfilter)
      olddat = ft_preproc_lowpassfilter(olddat, cfg.origfs, cfg.lpfreq, cfg.lpfiltord, cfg.lpfilttype);
    end

    % perform the resampling
    newtim = cfg.time{itr};
    if length(oldtim)>1
      newdat = interp1(oldtim', olddat', newtim', cfg.method, cfg.extrapval)';
    else
      newdat = repmat(olddat, [1 numel(newtim)]);
    end

    % add back the mean
    if ~strcmp(cfg.demean, 'yes')
      nsmp   = size(newdat, 2);
      newdat = newdat + bsl(:,ones(1,nsmp));
    end

    data.trial{itr} = newdat;
    data.time{itr}  = newtim;

  end % for itr
  ft_progress('close');

  % specify the new sampling frequency in the output
  t1 = cfg.time{1}(1);
  t2 = cfg.time{1}(2);
  data.fsample = 1/(t2-t1);

end % if usefsample or usetime

ft_info('original sampling rate = %d Hz\nnew sampling rate = %d Hz\n', cfg.origfs, data.fsample);

% convert back to input type if necessary
switch convert
  case 'timelock'
    data = ft_checkdata(data, 'datatype', 'timelock');
  otherwise
    % keep the output as it is
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   data
ft_postamble provenance data
ft_postamble history    data
ft_postamble savevar    data

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that decimates along the columns
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = my_decimate(x, varargin)
[n, m] = size(x);
% decimate the first column
y = decimate(x(:,1), varargin{:});
if m>1
  % increase the size of the output matrix
  y(:,m) = 0;
  % decimate the subsequent columns
  for i=2:m
    y(:,i) = decimate(x(:,i), varargin{:});
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that does a block-wise average along the columns
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = my_mean(x, r)
[n, m] = size(x);
n = n - mod(n,r);
x = x(1:n,:);
x = reshape(x, [r n/r m]);
y = mean(x, 1);
y = reshape(y, [n/r m]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that does a block-wise median along the columns
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = my_median(x, r)
[n, m] = size(x);
n = n - mod(n,r);
x = x(1:n,:);
x = reshape(x, [r n/r m]);
y = median(x, 1);
y = reshape(y, [n/r m]);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% This is a demonstration of Radon inversion from few angles using the Potts model.
% Tested on Mac OS 10.9, Matlab 2015a,
% If everything worked correctly a result as
% in results/plFigPottsRec7Angles.fig should appear
% Warning: This experiment can take more than one hour

% load Shepp-Logan
img = phantom(256);

% Radon transform with 7 equidistant projection angles
nAngles = 7;
theta = (1:nAngles)/nAngles * 180;
A = radonop(theta, size(img,1));
data = A * img;

% raw filtered backprojection and a regularized one
uFbp = iradon(data, theta, 'Linear', 'Ram-Lak', 1, size(img, 1));
uFbpReg = iradon(data, theta, 'Linear', 'Hamming', 0.3, size(img, 1));

% joint reconstruction and segmentation using the Potts model
gamma = 0.04;
uPotts = minL2iPotts2DADMM(data, gamma, A, 'verbose', 1, 'groundTruth', img, 'isotropic', 1 );

% show results
subplot(2,3,1)
imshow(img)
title('Original')
subplot(2,3,2)
imagesc(data)
title('Data (Sinogram)')
subplot(2,3,4)
imshow(uFbp)
title('FBP result')
subplot(2,3,5)
imshow(uFbpReg)
title('FBP with Hamming window')
subplot(2,3,6)
imshow(uPotts)
title('Potts result')

author: mstorath

repo: Pottslab

save_path: github-repos/MATLAB/mstorath-Pottslab


text:
function c=voigtUnMap(cVoigt)

if isvector(cVoigt) %assume that c is a 4th order tensor
    siz_c=[3 3];
    cVoigt(4:end)=(1/2).*cVoigt(4:end); %Undo doubling
else
    siz_c=[3 3 3 3];
end

[linearIndexVoigt,linearIndexFourthOrder]=tensor2voigtMap(zeros(siz_c));

switch class(cVoigt)
    case 'double'
        c=zeros(siz_c);
    case 'sym'
        c=sym(zeros(siz_c));
end

c(linearIndexFourthOrder)=cVoigt(linearIndexVoigt);

end

%%
% _*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 gegenbauer_rule ( order, alpha, a, b, filename )

%*****************************************************************************80
%
%% GEGENBAUER_RULE generates a Gauss-Gegenbauer rule.
%
%  Discussion:
%
%    This program computes a standard Gauss-Gegenbauer quadrature rule
%    and writes it to a file.
%
%    The user specifies:
%    * the ORDER (number of points) in the rule;
%    * the ALPHA parameter;
%    * A, the left endpoint;
%    * B, the right endpoint;
%    * FILENAME, the root name of the output files.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2010
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GEGENBAUER_RULE\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compute a Gauss-Gegenbauer rule for approximating\n' );
  fprintf ( 1, '    Integral ( A <= x <= B ) ((x-A)(B-X))^ALPHA f(x) dx\n' );
  fprintf ( 1, '  of order ORDER.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The user specifies ORDER, ALPHA, A, B, and FILENAME.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  ORDER is the number of points:\n' );
  fprintf ( 1, '  ALPHA is the exponent:\n' );
  fprintf ( 1, '  A is the left endpoint;\n' );
  fprintf ( 1, '  B is the right endpoint;\n' );
  fprintf ( 1, '  FILENAME is used to generate 3 files:\n' );
  fprintf ( 1, '    filename_w.txt - the weight file\n' );
  fprintf ( 1, '    filename_x.txt - the abscissa file.\n' );
  fprintf ( 1, '    filename_r.txt - the region file.\n' );
%
%  Initialize the parameters.
%
  beta = 0.0;
%
%  Get ORDER.
%
  if ( nargin < 1 )
    order = input ( '  Enter the rule order ORDER:  ' );
  elseif ( ischar ( order ) )
    order = str2num ( order );
  else
%
%  Get ALPHA.
%
  if ( nargin < 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  ALPHA is the exponent of ((x-A)(B-X)) in the integral:\n' );
    fprintf ( 1, '  Note that -1.0 < ALPHA is required.\n' );
    fprintf ( 1, '\n' );
    alpha = input ( '  Enter the value of ALPHA:  ' );
  elseif ( ischar ( alpha ) )
    alpha = str2num ( alpha );
  end
%
%  Get A.
%
  if ( nargin < 3 )
    a = input ( '  Enter the left endpoint A:  ' );
  elseif ( ischar ( a ) )
    a = str2num ( a );
  end
%
%  Get B.
%
  if ( nargin < 4 )
    b = input ( '  Enter the right endpoint B:  ' );
  elseif ( ischar ( b ) )
    b = str2num ( b );
  end
%
%  Get FILENAME.
%
  if ( nargin < 5 )
    fprintf ( 1,  '\n' );
    fprintf ( 1,  '  FILENAME is the ''root name'' of the quadrature files).\n' );
    filename = input ( '  Enter FILENAME as a quoted string:  ' );
  end
%
%  Input summary.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  ORDER = %d\n', order );
  fprintf ( 1, '  ALPHA = %f\n', alpha );
  fprintf ( 1, '  A = %f\n', a );
  fprintf ( 1, '  B = %f\n', b );
  fprintf ( 1, '  FILENAME = "%s".\n', filename );
%
%  Construct the rule.
%
  kind = 3;
  [ x, w ] = cgqf ( order, kind, alpha, beta, a, b );
%
%  Write the rule.
%
  r = [ a, b ]';
  rule_write ( order, filename, x, w, r );
%
%  Terminate.
%
  fprintf ( 1,  '\n' );
  fprintf ( 1,  'GEGENBAUER_RULE:\n' );
  fprintf ( 1,  '  Normal end of execution.\n' );
  fprintf ( 1,  '\n' );
  timestamp ( );

  return
end
function [ t, wts ] = cdgqf ( nt, kind, alpha, beta )

%*****************************************************************************80
%
%% CDGQF computes a Gauss quadrature formula with default A, B and simple knots.
%
%  Discussion:
%
%    This routine computes all the knots and weights of a Gauss quadrature
%    formula with a classical weight function with default values for A and B,
%    and only simple knots.
%
%    There are no moments checks and no printing is done.
%
%    Use routine EIQFS to evaluate a quadrature computed by CGQFS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev,            (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,inf)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-inf,inf)  |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,inf)     (x-a)^alpha*(x+b)^beta
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Output, real T(NT), the knots.
%
%    Output, real WTS(NT), the weights.
%
  parchk ( kind, 2 * nt, alpha, beta );
%
%  Get the Jacobi matrix and zero-th moment.
%
  [ aj, bj, zemu ] = class_matrix ( kind, nt, alpha, beta );
%
%  Compute the knots and weights.
%
  [ t, wts ] = sgqf ( nt, aj, bj, zemu );

  return
end
function [ t, wts ] = cgqf ( nt, kind, alpha, beta, a, b )

%*****************************************************************************80
%
%% CGQF computes knots and weights of a Gauss quadrature formula.
%
%  Discussion:
%
%    The user may specify the interval (A,B).
%
%    Only simple knots are produced.
%
%    The user may request that the routine print the knots and weights,
%    and perform a moment check.
%
%    Use routine EIQFS to evaluate this quadrature formula.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev Type 1,     (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,+oo)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-oo,+oo)   |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,+oo)     (x-a)^alpha*(x+b)^beta
%    9, Chebyshev Type 2,     (a,b)       ((b-x)*(x-a))^(+0.5)
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Input, real A, B, the interval endpoints.
%
%    Output, real T(NT), the knots.
%
%    Output, real WTS(NT), the weights.
%

%
%  Compute the Gauss quadrature formula for default values of A and B.
%
  [ t, wts ] = cdgqf ( nt, kind, alpha, beta );
%
%  All knots have multiplicity = 1.
%
  mlt = zeros(nt,1);
  mlt(1:nt) = 1;
%
%  NDX(I) = I.
%
  ndx = ( 1 : nt );
%
%  Scale the quadrature rule.
%
  [ t, wts ] = scqf ( nt, t, mlt, wts, nt, ndx, kind, alpha, beta, a, b );

  return
end
function [ aj, bj, zemu ] = class_matrix ( kind, m, alpha, beta )

%*****************************************************************************80
%
%% CLASS_MATRIX computes the Jacobi matrix for a quadrature rule.
%
%  Discussion:
%
%    This routine computes the diagonal AJ and subdiagonal BJ
%    elements of the order M tridiagonal symmetric Jacobi matrix
%    associated with the polynomials orthogonal with respect to
%    the weight function specified by KIND.
%
%    For weight functions 1-7, M elements are defined in BJ even
%    though only M-1 are needed.  For weight function 8, BJ(M) is
%    set to zero.
%
%    The zero-th moment of the weight function is returned in ZEMU.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev,            (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,inf)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-inf,inf)  |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,inf)     (x-a)^alpha*(x+b)^beta
%
%    Input, integer M, the order of the Jacobi matrix.
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Output, real AJ(M), BJ(M), the diagonal and subdiagonal
%    of the Jacobi matrix.
%
%    Output, real ZEMU, the zero-th moment.
%
  temp = eps;

  parchk ( kind, 2 * m - 1, alpha, beta );

  temp2 = 0.5;

  if ( 500.0 * temp < abs ( ( gamma ( temp2 ) )^2 - pi ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CLASS - Fatal error!\n' );
    fprintf ( 1, '  Gamma function does not match machine parameters.\n' );
    error ( 'CLASS - Fatal error!' );
  end

  bj = zeros(m,1);
  aj = zeros(m,1);

  if ( kind == 1 )

    ab = 0.0;

    zemu = 2.0 / ( ab + 1.0 );

    aj(1:m) = 0.0;

    for i = 1 : m
      abi = i + ab * mod ( i, 2 );
      abj = 2 * i + ab;
      bj(i) = abi * abi / ( abj * abj - 1.0 );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 2 )

    zemu = pi;

    aj(1:m) = 0.0;

    bj(1) =  sqrt ( 0.5 );
    bj(2:m) = 0.5;

  elseif ( kind == 3 )

    ab = alpha * 2.0;
    zemu = 2.0^( ab + 1.0 ) * gamma ( alpha + 1.0 )^2 ...
      / gamma ( ab + 2.0 );

    aj(1:m) = 0.0;
    bj(1) = 1.0 / ( 2.0 * alpha + 3.0 );
    for i = 2 : m
      bj(i) = i * ( i + ab ) / ( 4.0 * ( i + alpha )^2 - 1.0 );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 4 )

    ab = alpha + beta;
    abi = 2.0 + ab;
    zemu = 2.0^( ab + 1.0 ) * gamma ( alpha + 1.0 ) ...
      * gamma ( beta + 1.0 ) / gamma ( abi );
    aj(1) = ( beta - alpha ) / abi;
    bj(1) = 4.0 * ( 1.0 + alpha ) * ( 1.0 + beta ) ...
      / ( ( abi + 1.0 ) * abi * abi );
    a2b2 = beta * beta - alpha * alpha;

    for i = 2 : m
      abi = 2.0 * i + ab;
      aj(i) = a2b2 / ( ( abi - 2.0 ) * abi );
      abi = abi^2;
      bj(i) = 4.0 * i * ( i + alpha ) * ( i + beta ) * ( i + ab ) ...
        / ( ( abi - 1.0 ) * abi );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 5 )

    zemu = gamma ( alpha + 1.0 );

    for i = 1 : m
      aj(i) = 2.0 * i - 1.0 + alpha;
      bj(i) = i * ( i + alpha );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 6 )

    zemu = gamma ( ( alpha + 1.0 ) / 2.0 );

    aj(1:m) = 0.0;

    for i = 1 : m
      bj(i) = ( i + alpha * mod ( i, 2 ) ) / 2.0;
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 7 )

    ab = alpha;
    zemu = 2.0 / ( ab + 1.0 );

    aj(1:m) = 0.0;

    for i = 1 : m
      abi = i + ab * mod(i,2);
      abj = 2 * i + ab;
      bj(i) = abi * abi / ( abj * abj - 1.0 );
    end
    bj(1:m) =  sqrt ( bj(1:m) );

  elseif ( kind == 8 )

    ab = alpha + beta;
    zemu = gamma ( alpha + 1.0 ) * gamma ( - ( ab + 1.0 ) ) ...
      / gamma ( - beta );
    apone = alpha + 1.0;
    aba = ab * apone;
    aj(1) = - apone / ( ab + 2.0 );
    bj(1) = - aj(1) * ( beta + 1.0 ) / ( ab + 2.0 ) / ( ab + 3.0 );
    for i = 2 : m
      abti = ab + 2.0 * i;
      aj(i) = aba + 2.0 * ( ab + i ) * ( i - 1 );
      aj(i) = - aj(i) / abti / ( abti - 2.0 );
    end

    for i = 2 : m - 1
      abti = ab + 2.0 * i;
      bj(i) = i * ( alpha + i ) / ( abti - 1.0 ) * ( beta + i ) ...
        / ( abti^2 ) * ( ab + i ) / ( abti + 1.0 );
    end

    bj(m) = 0.0;
    bj(1:m) =  sqrt ( bj(1:m) );

  end

  return
end
function [ d, z ] = imtqlx ( n, d, e, z )

%*****************************************************************************80
%
%% IMTQLX diagonalizes a symmetric tridiagonal matrix.
%
%  Discussion:
%
%    This routine is a slightly modified version of the EISPACK routine to
%    perform the implicit QL algorithm on a symmetric tridiagonal matrix.
%
%    The authors thank the authors of EISPACK for permission to use this
%    routine.
%
%    It has been modified to produce the product Q' * Z, where Z is an input
%    vector and Q is the orthogonal matrix diagonalizing the input matrix.
%    The changes consist (essentialy) of applying the orthogonal transformations
%    directly to Z as they are generated.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%    Roger Martin, James Wilkinson,
%    The Implicit QL Algorithm,
%    Numerische Mathematik,
%    Volume 12, Number 5, December 1968, pages 377-383.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real D(N), the diagonal entries of the matrix.
%
%    Input, real E(N), the subdiagonal entries of the
%    matrix, in entries E(1) through E(N-1).
%
%    Input, real Z(N), a vector to be operated on.
%
%    Output, real D(N), the diagonal entries of the diagonalized matrix.
%
%    Output, real Z(N), the value of Q' * Z, where Q is the matrix that
%    diagonalizes the input symmetric tridiagonal matrix.
%
  itn = 30;

  prec = eps;

  if ( n == 1 )
    return
  end

  e(n) = 0.0;

  for l = 1 : n

    j = 0;

    while ( 1 )

      for m = l : n

        if ( m == n )
          break
        end

        if ( abs ( e(m) ) <= prec * ( abs ( d(m) ) + abs ( d(m+1) ) ) )
          break
        end

      end

      p = d(l);

      if ( m == l )
        break
      end

      if ( j == itn )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'IMTQLX - Fatal error!\n' );
        fprintf ( 1, '  Iteration limit exceeded.\n' );
        error ( 'IMTQLX - Fatal error!' );
      end

      j = j + 1;
      g = ( d(l+1) - p ) / ( 2.0 * e(l) );
      r =  sqrt ( g * g + 1.0 );
      g = d(m) - p + e(l) / ( g + r8_sign ( g ) * abs ( r ) );
      s = 1.0;
      c = 1.0;
      p = 0.0;
      mml = m - l;

      for ii = 1 : mml

        i = m - ii;
        f = s * e(i);
        b = c * e(i);

        if ( abs ( f ) >= abs ( g ) )
          c = g / f;
          r =  sqrt ( c * c + 1.0 );
          e(i+1) = f * r;
          s = 1.0 / r;
          c = c * s;
        else
          s = f / g;
          r =  sqrt ( s * s + 1.0 );
          e(i+1) = g * r;
          c = 1.0 / r;
          s = s * c;
        end

        g = d(i+1) - p;
        r = ( d(i) - g ) * s + 2.0 * c * b;
        p = s * r;
        d(i+1) = g + p;
        g = c * r - b;
        f = z(i+1);
        z(i+1) = s * z(i) + c * f;
        z(i) = c * z(i) - s * f;

      end

      d(l) = d(l) - p;
      e(l) = g;
      e(m) = 0.0;

    end

  end

  for ii = 2 : n

     i = ii - 1;
     k = i;
     p = d(i);

     for j = ii : n
       if ( d(j) < p )
         k = j;
         p = d(j);
       end
     end

     if ( k ~= i )
       d(k) = d(i);
       d(i) = p;
       p = z(i);
       z(i) = z(k);
       z(k) = p;
     end

  end

  return
end
function parchk ( kind, m, alpha, beta )

%*****************************************************************************80
%
%% PARCHK checks parameters ALPHA and BETA for classical weight functions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 January 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev,            (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,inf)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-inf,inf)  |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,inf)     (x-a)^alpha*(x+b)^beta
%
%    Input, integer M, the order of the highest moment to
%    be calculated.  This value is only needed when KIND = 8.
%
%    Input, real ALPHA, BETA, the parameters, if required
%    by the value of KIND.
%
  if ( kind <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PARCHK - Fatal error!\n' );
    fprintf ( 1, '  KIND <= 0.\n' );
    error ( 'PARCHK - Fatal error!' );
  end
%
%  Check ALPHA for Gegenbauer, Jacobi, Laguerre, Hermite, Exponential.
%
  if ( 3 <= kind && alpha <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PARCHK - Fatal error!\n' );
    fprintf ( 1, '  3 <= KIND and ALPHA <= -1.\n' );
    error ( 'PARCHK - Fatal error!' );
  end
%
%  Check BETA for Jacobi.
%
  if ( kind == 4 && beta <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PARCHK - Fatal error!\n' );
    fprintf ( 1, '  KIND == 4 and BETA <= -1.0.\n' );
    error ( 'PARCHK - Fatal error!' );
  end
%
%  Check ALPHA and BETA for rational.
%
  if ( kind == 8 )
    tmp = alpha + beta + m + 1.0;
    if ( 0.0 <= tmp || tmp <= beta )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PARCHK - Fatal error!\n' );
      fprintf ( 1, '  KIND == 8 but condition on ALPHA and BETA fails.\n' );
      error ( 'PARCHK - Fatal error!' );
    end
  end

  return
end
function value = r8_sign ( x )

%*****************************************************************************80
%
%% R8_SIGN returns the sign of an R8.
%
%  Discussion:
%
%    The value is +1 if the number is positive or zero, and it is -1 otherwise.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the number whose sign is desired.
%
%    Output, real VALUE, the sign of X.
%
  if ( 0 <= x )
    value = +1.0;
  else
    value = -1.0;
  end

  return
end
function r8mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% R8MAT_WRITE writes an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE - Error!' );
  end
%
%  Write the data.
%
%  For smaller data files, and less precision, try:
%
%     fprintf ( output_unit, '  %14.6f', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %24.16f', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function rule_write ( order, filename, x, w, r )

%*****************************************************************************80
%
%% RULE_WRITE writes a quadrature rule to a file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ORDER, the order of the rule.
%
%    Input, string FILENAME, specifies the output files.
%    write files 'filename_w.txt', 'filename_x.txt', 'filename_r.txt' defining
%    weights, abscissas, and region.
%
%    Input, real X(ORDER), the abscissas.
%
%    Input, real W(ORDER), the weights.
%
%    Input, real R(2), the region.
%
  filename_x = strcat ( filename, '_x.txt' );
  filename_w = strcat ( filename, '_w.txt' );
  filename_r = strcat ( filename, '_r.txt' );

  fprintf ( 1, '\n' );
  fprintf ( 1,'  Creating quadrature files.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  "Root" file name is   "%s".\n', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Weight file will be   "%s".\n', filename_w );
  fprintf ( 1, '  Abscissa file will be "%s".\n', filename_x );
  fprintf ( 1, '  Region file will be   "%s".\n', filename_r );

  r8mat_write ( filename_w, 1, order, w' );
  r8mat_write ( filename_x, 1, order, x' );
  r8mat_write ( filename_r, 1, 2,     r' );

  return
end
function [ t, wts ] = scqf ( nt, t, mlt, wts, nwts, ndx, kind, alpha, ...
  beta, a, b )

%*****************************************************************************80
%
%% SCQF scales a quadrature formula to a nonstandard interval.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, real T(NT), the original knots.
%
%    Input, integer MLT(NT), the multiplicity of the knots.
%
%    Input, real WTS(NWTS), the weights.
%
%    Input, integer NWTS, the number of weights.
%
%    Input, integer NDX(NT), used to index the array WTS.
%    For more details see the comments in CAWIQ.
%
%    Input, integer KIND, the rule.
%    1, Legendre,             (a,b)       1.0
%    2, Chebyshev Type 1,     (a,b)       ((b-x)*(x-a))^(-0.5)
%    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
%    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
%    5, Generalized Laguerre, (a,+oo)     (x-a)^alpha*exp(-b*(x-a))
%    6, Generalized Hermite,  (-oo,+oo)   |x-a|^alpha*exp(-b*(x-a)^2)
%    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
%    8, Rational,             (a,+oo)     (x-a)^alpha*(x+b)^beta
%    9, Chebyshev Type 2,     (a,b)       ((b-x)*(x-a))^(+0.5)
%
%    Input, real ALPHA, the value of Alpha, if needed.
%
%    Input, real BETA, the value of Beta, if needed.
%
%    Input, real A, B, the interval endpoints.
%
%    Output, real T(NT), the scaled knots.
%
%    Output, real WTS(NWTS), the scaled weights.
%
  temp = eps;

  parchk ( kind, 1, alpha, beta )

  if ( kind == 1 )

    al = 0.0;
    be = 0.0;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 2 )

    al = -0.5;
    be = -0.5;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 3 )

    al = alpha;
    be = alpha;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 4 )

    al = alpha;
    be = beta;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 5 )

    if ( b <= 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  B <= 0.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = a;
    slp = 1.0 / b;
    al = alpha;
    be = 0.0;

  elseif ( kind == 6 )

    if ( b <= 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  B <= 0.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = a;
    slp = 1.0 / sqrt ( b );
    al = alpha;
    be = 0.0;

  elseif ( kind == 7 )

    al = alpha;
    be = 0.0;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  elseif ( kind == 8 )

    if ( a + b <= 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  A + B <= 0.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = a;
    slp = a + b;
    al = alpha;
    be = beta;

  elseif ( kind == 9 )

    al = 0.5;
    be = 0.5;

    if ( abs ( b - a ) <= temp )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SCQF - Fatal error!\n' );
      fprintf ( 1, '  |B - A| too small.\n' );
      fprintf ( 1, '  A = %f\n', a );
      fprintf ( 1, '  B = %f\n', b );
      error ( 'SCQF - Fatal error!' );
    end

    shft = ( a + b ) / 2.0;
    slp = ( b - a ) / 2.0;

  end

  p = slp^( al + be + 1.0 );

  for k = 1 : nt

    t(k) = shft + slp * t(k);
    l = abs ( ndx(k) );

    if ( l ~= 0 )
      tmp = p;
      for i = l : l + mlt(k) - 1
        wts(i) = wts(i) * tmp;
        tmp = tmp * slp;
      end
    end

  end

  return
end
function [ t, wts ] = sgqf ( nt, aj, bj, zemu )

%*****************************************************************************80
%
%% SGQF computes knots and weights of a Gauss Quadrature formula.
%
%  Discussion:
%
%    This routine computes all the knots and weights of a Gauss quadrature
%    formula with simple knots from the Jacobi matrix and the zero-th
%    moment of the weight function, using the Golub-Welsch technique.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 February 2010
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer NT, the number of knots.
%
%    Input, real AJ(NT), the diagonal of the Jacobi matrix.
%
%    Input, real BJ(NT), the subdiagonal of the Jacobi
%    matrix, in entries 1 through NT-1.  On output, BJ has been overwritten.
%
%    Input, real ZEMU, the zero-th moment of the weight function.
%
%    Output, real T(NT), the knots.
%
%    Output, real WTS(NT), the weights.
%

%
%  Exit if the zero-th moment is not positive.
%
  if ( zemu <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SGQF - Fatal error!\n' );
    fprintf ( 1, '  ZEMU <= 0.\n' );
    error ( 'SGQF - Fatal error!' );
  end
%
%  Set up vectors for IMTQLX.
%
  wts = zeros ( nt, 1 );

  wts(1) = sqrt ( zemu );
  wts(2:nt) = 0.0;
%
%  Diagonalize the Jacobi matrix.
%
  [ t, wts ] = imtqlx ( nt, aj, bj, wts );

  wts(1:nt) = wts(1:nt).^2;

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [x, y, z] = torus (a, n, r, kpi)
% TORUS Generate a torus.
% torus (r, n, a, kpi) generates a plot of a
% torus with central radius a and
% lateral radius r.
% n controls the number of facets
% on the surface.
% kpi makes it possible to draw a whole torus,
% or e.g. half of it.
%
% This script is a modification of a
% program from:

% MATLAB Primer, 6th Edition
% Kermit Sigmon and Timothy A. Davis
% Section 11.5, page 65.

theta = -pi * (0:2:kpi*n)/n ;
phi = 2*pi* (0:2:n)'/n ;
x = (a + r*cos(phi)) * cos(theta) ;
y = (a + r*cos(phi)) * sin(theta) ;
z = r * sin(phi) * ones(size(theta)) ;

author: Sable

repo: mcbench-benchmarks


text:
function [xi,w]=seventhOrder2DCubPoints()
%%SEVENTHORDER2DCUBPOINTS Generate seventh-order cubature points
%               for integration over a 2-dimensional cube with bounds of
%               (-1,-1), (-1,1), (1,1), and (1,-1).
%
%INPUTS: None
%
%OUTPUTS: xi This is a 2XnumCubPoints set of points for the standard
%            square.
%          w A 1XnumCubPoints set of cubature weights. This sums to the
%            volume of the standard square (4).
%
%This function implements the points given in [1] (12 points).
%
%EXAMPLE:
%We compare a 6th-order moment computed using these cubature points
%to one computed using monomialIntCube (a 7th order moment would have just
%been 0). The results are the same within typical finite precision limits.
% [xi,w]=seventhOrder2DCubPoints();
% alpha=[2;4];
% theMoment=findMomentFromSamp(alpha,xi,w);
% intVal=monomialIntCube(alpha);
% RelErr=(theMoment-intVal)/intVal
%
%REFERENCES:
%[1] F. D. Witherden and P. E. Vincent, "On the identification of symmetric
%    quadrature rules for finite element methods," Computer and Mathematics
%    with Applications, vol. 69, no. 10, pp. 1232-1241, May 2015.
%
%October 2022 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

M=[0.92582009977255146156656677658399952253,                                          0,   0.24197530864197530864197530864197530864;
                                        0,   0.92582009977255146156656677658399952253,   0.24197530864197530864197530864197530864;
-0.92582009977255146156656677658399952253,                                          0,   0.24197530864197530864197530864197530864;
                                        0,  -0.92582009977255146156656677658399952253,   0.24197530864197530864197530864197530864;
  0.8059797829185987437078561813507442463,    0.8059797829185987437078561813507442463,   0.23743177469063023421810525931129352533;
  0.8059797829185987437078561813507442463,   -0.8059797829185987437078561813507442463,   0.23743177469063023421810525931129352533;
 -0.8059797829185987437078561813507442463,    0.8059797829185987437078561813507442463,   0.23743177469063023421810525931129352533;
 -0.8059797829185987437078561813507442463,   -0.8059797829185987437078561813507442463,   0.23743177469063023421810525931129352533;
  0.3805544332083156563791063590863941355,    0.3805544332083156563791063590863941355,   0.52059291666739445713991943204673116603;
  0.3805544332083156563791063590863941355,   -0.3805544332083156563791063590863941355,   0.52059291666739445713991943204673116603;
 -0.3805544332083156563791063590863941355,    0.3805544332083156563791063590863941355,   0.52059291666739445713991943204673116603;
 -0.3805544332083156563791063590863941355,   -0.3805544332083156563791063590863941355,   0.52059291666739445713991943204673116603];
w=M(:,3);
xi=M(:,1:2)';

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:
%IMM_UPDATE  Interacting Multiple Model (IMM) Filter update step
%
% Syntax:
%   [X_i,P_i,MU,X,P] = IMM_UPDATE(X_p,P_p,c_j,ind,dims,Y,H,R)
%
% In:
%   X_p  - Cell array containing N^j x 1 mean state estimate vector for
%          each model j after prediction step
%   P_p  - Cell array containing N^j x N^j state covariance matrix for
%          each model j after prediction step
%   c_j  - Normalizing factors for mixing probabilities
%   ind  - Indices of state components for each model as a cell array
%   dims - Total number of different state components in the combined system
%   Y    - Dx1 measurement vector.
%   H    - Measurement matrices for each model as a cell array.
%   R    - Measurement noise covariances for each model as a cell array.
%
% Out:
%   X_i  - Updated state mean estimate for each model as a cell array
%   P_i  - Updated state covariance estimate for each model as a cell array
%   MU   - Estimated probabilities of each model
%   X    - Combined state mean estimate
%   P    - Combined state covariance estimate
%
% Description:
%   IMM filter measurement update step.
%
% See also:
%   IMM_PREDICT, IMM_SMOOTH, IMM_FILTER

% History:
%   01.11.2007 JH The first official version.
%
% Copyright (C) 2007 Jouni Hartikainen
%
% $Id: imm_update.m 111 2007-11-01 12:09:23Z jmjharti $
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.

function [X_i,P_i,MU,X,P] = imm_update(X_p,P_p,c_j,ind,dims,Y,H,R)
    % Number of models
    m = length(X_p);

    % Space for update state mean, covariance and likelihood of measurements
    X_i = cell(1,m);
    P_i = cell(1,m);
    lambda = zeros(1,m);

    % Update for each model
    for i = 1:m
        % Update the state estimates
        [X_i{i}, P_i{i}, K, IM, IS, lambda(i)] = kf_update(X_p{i},P_p{i},Y,H{i},R{i});
    end

    % Calculate the model probabilities
    MU = zeros(1,m);
    c = sum(lambda.*c_j);
    MU = c_j.*lambda/c;

    % Output the combined updated state mean and covariance, if wanted.
    if nargout > 3
        % Space for estimates
        X = zeros(dims,1);
        P = zeros(dims,dims);
        % Updated state mean
        for i = 1:m
            X(ind{i}) = X(ind{i}) + MU(i)*X_i{i};
        end
        % Updated state covariance
        for i = 1:m
            P(ind{i},ind{i}) = P(ind{i},ind{i}) + MU(i)*(P_i{i} + (X_i{i}-X(ind{i}))*(X_i{i}-X(ind{i}))');
        end
    end

author: EEA-sensors

repo: ekfukf

save_path: github-repos/MATLAB/EEA-sensors-ekfukf


text:
clear all;
close all;
clc;

image = 'cow';
% image = 'bungee';
% image = 'man';

Aorg = imread([image '.png']);
Morg = imread([image '-mask.png']);

tic
A=patch_inpaint(Aorg,Morg);
toc

imshow(A)


text:
function ADMINrenamefilenames(HOMEANNOTATIONS, Folder)
%
% Writes in the filename field of the xml file the nane of the filename.
%
% This script is useful when changing image and annotation file names. This function guaratees consistency.

% Create list of folders
if nargin == 1
    Folder = folderlist(HOMEANNOTATIONS);
    Folder = [{''} Folder];
end

Nfolders = length(Folder);

for n = 1:Nfolders
    f = Folder{n};

    if ~isempty(fullfile(HOMEANNOTATIONS, f))
        files = dir(fullfile(HOMEANNOTATIONS, f, '*.xml'));

        Nfiles = length(files);
        for i = 1:Nfiles
            newname = files(i).name;
            newname = strrep(newname, '.xml', '.jpg');

            % load annotation
            dest = fullfile(HOMEANNOTATIONS, f, files(i).name);
            xml = loadxml(dest);

            % change inside annotation file
            ii = strfind(xml, '<filename>');
            jj = strfind(xml, '</filename>');
            oldname = strtrim(xml(ii+10:jj-1));

            if ~strcmp(oldname, newname)
                disp('rename filename field inside xml')

                xml = [xml(1:ii+9) newname xml(jj:end)];

                % save annotation
                fid = fopen(dest, 'w');
                fprintf(fid, xml');
                fclose(fid);
            end
        end
    end

end

function xml = loadxml(filename)
[fid, message] = fopen(filename,'r');
if fid == -1; error(message); end
xml = fread(fid, 'uint8=>char');
fclose(fid);
xml = xml';

author: CSAILVision

repo: LabelMeToolbox


text:
classdef HybridSystem < handle & matlab.mixin.Copyable
    % HybridSystem defines a hybrid dynamical system that has both
    % continuous and discrete dynamics, such as bipedal locomotion.
    %
    % This class provides basic elements and functionalities of a hybrid
    % dynamicsl system. The mathematical definition of the hybrid system is
    % given as
    % \f{eqnarray*}{
    % \mathscr{HC} = \{\Gamma, \mathcal{D}, U, S, \Delta, FG\}
    % \f}
    %
    % The implementation of the hybrid system class is heavily based on the
    % Matlab's digraph data type, with wrapper functions with additional
    % validation.
    %
    % @author Ayonga Hereid @date 2016-09-26
    %
    % Copyright (c) 2016, AMBER Lab
    % All right reserved.
    %
    % Redistribution and use in source and binary forms, with or without
    % modification, are permitted only in compliance with the BSD 3-Clause
    % license, see
    % http://www.opensource.org/licenses/bsd-license.php

    %% public properties
    properties (GetAccess = public, SetAccess = protected)
        % This is the name of the object that gives the object an universal
        % identification
        %
        % @type char @default ''
        Name

        % The directed graph that describes the hybrid dynamical system
        % structure
        %
        % @type digraph
        Gamma

        % The option of the hybrid system model
        %
        % @type struct
        Options

    end

    properties (Dependent)
        % A structure describes the name, type, attributes, and
        % default value of each properties of the vertex
        %
        % @type struct
        VertexProperties

        % A structure describes the name, type, attributes, and
        % default value of each properties of the edge
        %
        % @type struct
        EdgeProperties
    end
    methods
        function VertexProperties = get.VertexProperties(~)

            VertexProperties = struct();
            VertexProperties.Name =  {'Domain','Control','Param'};
            VertexProperties.Type = {{'ContinuousDynamics'},{'Controller'},{'struct'}};
            VertexProperties.Attribute = {{},{},{}};
            VertexProperties.DefaultValue =  {{[]},{[]},{[]}};
        end

        function EdgeProperties = get.EdgeProperties(~)

            EdgeProperties = struct();
            EdgeProperties.Name =  {'Guard', 'Param', 'Weights'};
            EdgeProperties.Type = {{'DiscreteDynamics'}, {'struct'}, {'numeric'}};
            EdgeProperties.Attribute = {{}, {}, {'scalar'}};
            EdgeProperties.DefaultValue = {{[]},{[]}, NaN};
        end
    end

    %% Public methods
    methods (Access = public)
        function obj = HybridSystem(name, varargin)
            % the default calss constructor
            %
            % Parameters:
            % name: the name of the hybrid system model @type char
            % varargin: variable parameter input argument @type varargin
            %
            % Return values:
            % obj: the class object

            assert(ischar(name), 'The object name must be a character vector.');
            obj.Name = name;

            % initialize an empty directed graph with specified properties
            obj.Gamma = digraph;
            % add a dummy node and remove it to initialize the properties
            obj = addVertex(obj,'dummy');
            obj = rmVertex(obj,'dummy');

            % load default options
            obj.Options = struct;
            obj.Options.Logger = 'SimLogger';
            obj.Options.OdeSolver = @ode45;

            % parse the input options
            if nargin > 1
                % create a structure from the input argument
                opts = struct(varargin{:});
                % overwrite the default options
                obj.setOption(opts);
            end

        end

        function obj = setOption(obj, varargin)
            % Sets the object options, and return the complete list of option
            % structure including unchanged default options.
            %
            % Parameters:
            % varargin: name-value pairs of option field value
            %
            % Return values:
            % options: the complete list of final option structure

            new_opts = struct(varargin{:});

            obj.Options = struct_overlay(obj.Options, new_opts,{'AllowNew',true});

        end

        function ret = isdag(obj)
            % returns true if the directed graph is a directed acyclic
            % graph
            %
            % @see ISDAG

            ret = isdag(obj.Gamma);

        end
        function ret = isDirectedCycle(obj)
            % returns true if the underlying directed graph is a simple
            % directed cycle.
            %
            % A simple directed cycle is a directed graph that has uniform
            % in-degree 1 and uniform out-degree 1.

            g = obj.Gamma;

            ret = all(indegree(g)==1) && all(outdegree(g)==1);

        end

        function sys = subGraph(obj, nodeIDs)
            % extract the subgraph of the hybrid system to create a new
            % hybrid system object with the same name
            %
            % Parameters:
            % nodeIDs: the node ids of the subgraph
            %
            % Return values:
            % sys: a new HybridSystem object with the subgraph specified by
            % the nodeIDs @type HybridSystem

            sys = HybridSystem(obj.Name);

            sub_g = subgraph(obj.Gamma, nodeIDs);

            sys.Gamma = sub_g;
        end

    end

    %% methods defined in separate files
    methods
        obj = addVertex(obj, varargin);

        obj = addEdge(obj, varargin);

        obj = rmVertex(obj, vertex_name);

        obj = rmEdge(obj, s, t);

        obj = setEdgeProperties(obj, s, t, varargin);

        obj = setVertexProperties(obj, vertex, varargin);

        obj = simulate(obj, t0, x0, tf, logger, options, varargin);

        obj = compile(obj, export_path, varargin);

        obj = saveExpression(obj, export_path, varargin);
    end

    %% Private methods
    methods (Static, Access=private)
        function validatePropAttribute(name, value, type, attribute)

            if iscell(value)
                cellfun(@(x)validateattributes(x,type,attribute,...
                    'HybridSystem', name), value);
            elseif isnumeric(value)
                % if a property is a numeric value, it must be a row vector
                for i=1:size(value,1)
                    validateattributes(value(i,:),type,attribute,'HybridSystem', name);
                end
            else
                error('The input argument must be a cell array of objects.');
            end

        end
    end

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function n2 = sp_dist2(x, c)
% DIST2	Calculates squared distance between two sets of points.
% Adapted from Netlab neural network software:
% http://www.ncrg.aston.ac.uk/netlab/index.php
%
%	Description
%	D = DIST2(X, C) takes two matrices of vectors and calculates the
%	squared Euclidean distance between them.  Both matrices must be of
%	the same column dimension.  If X has M rows and N columns, and C has
%	L rows and N columns, then the result has M rows and L columns.  The
%	I, Jth entry is the  squared distance from the Ith row of X to the
%	Jth row of C.
%
%
%	Copyright (c) Ian T Nabney (1996-2001)

[ndata, dimx] = size(x);
[ncentres, dimc] = size(c);
if dimx ~= dimc
	error('Data dimension does not match dimension of centres')
end

n2 = (ones(ncentres, 1) * sum((x.^2)', 1))' + ...
  ones(ndata, 1) * sum((c.^2)',1) - ...
  2.*(x*(c'));

% Rounding errors occasionally cause negative entries in n2
if any(any(n2<0))
  n2(n2<0) = 0;
end

author: adikhosla

repo: feature-extraction


text:
function D = sampling_func(f, domain, gridsize)
%SAMPLING_FUNC Sampling a (multivariate) function
%	D = SAMPLING_FUNC(f, domain, grid)
%
%	f      - function handle with a real vector argument with N elements
%	domain - [min1 max1;... minN maxN] intervals for each element
%	grid   - number of sampling grid points for each element
%
%	D      - N-dimensional array of the sampled data
%
%   eg.:   sampling_func(@(x) x(1)+x(2), [-1 1; 0 3], [7 5])

N = size(domain, 1);

if length(gridsize) == 1
	gridsize = gridsize * ones(N,1);
else
	if not(length(gridsize) == N)
		error 'length of "gridsize" and "domain" must be the same.'
	end
	gridsize = gridsize(:);
end

if N > 0
	% allocation of sample array
	siz = prod(gridsize);
	D = zeros(siz,1);
	a = domain(:,1);
	b = domain(:,2);
	step = (b - a) ./ (gridsize - 1);

	% sampling
	z = ones(N,1);
	z(1) = 0;
	for k = 1:siz
		% next gridsize point index
		z = nexti(gridsize, z);

		% arguments
		p = a + step .* (z - 1);

		% sampling the model at the gridsize point
		D(k) = f(p);
	end

	if N > 1
		% reshape to proper size
		D = reshape(D, gridsize');
	end
else
	D = f([]);
end

% next gridsize point index
function n = nexti(gridsize, n)
i = 1;
n(i) = n(i) + 1;
while n(i) > gridsize(i)
	n(i) = 1;
	i = i + 1;
	n(i) = n(i) + 1;
end


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: September 9th, 2016
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting %%	lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us %to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the geometry and prints associated input files
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Make_Scallop_Geo_and_Input_Files()

%
% Grid Parameters (MAKE SURE MATCHES IN input2d !!!)
%
Nx =  256;        % # of Eulerian Grid Pts. in x-Direction (MUST BE EVEN!!!)
Lx = 1.0;        % Length of Eulerian Grid in x-Direction
dx = Lx/Nx;      % Grid spatial resolution

%
% Immersed Structure Geometric / Dynamic Parameters %
%
ds = 0.5*dx;              % Lagrangian Pt. Spacing (2x resolution of Eulerian grid)
struct_name = 'scallop'; % Name for .vertex, .spring, etc files. (must match what's in 'input2d')

% Call function to construct geometry
[xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(ds,Lx);

half_len = (length(xLag)-1)/2;   % Computes # of Lag Pts on each "Arm"
ind_off = ceil( 0.4*half_len ); % Computes # of lags from center pt for muscle placement
b_ind = half_len - ind_off;      % Bottom index
t_ind = half_len+2 + ind_off;    % Top index

%b_ind = 1;
%t_ind = length(xLag);
m_dist = sqrt( ( xLag(b_ind)-xLag(t_ind) )^2 + ( yLag(b_ind)-yLag(t_ind) )^2   );

% Test Muscle Placement
plot( xLag( b_ind ), yLag( b_ind ), 'm*'); hold on;
plot( xLag( t_ind ), yLag( t_ind ), 'm*'); hold on;

% Prints .vertex file!
print_Lagrangian_Vertices(xLag,yLag,struct_name);

% Prints .spring file!
k_Spring = 2e6;                    % Spring stiffness (does not need to be equal for all springs)
ds_Rest = ds;                        % Spring resting length (does not need to be equal for all springs)
print_Lagrangian_Springs(xLag,yLag,k_Spring,ds_Rest,struct_name);

% Prints .beam file!
k_Beam = 1e11;                      % Beam Stiffness (does not need to be equal for all beams)
C = compute_Curvatures(xLag,yLag); % Computes curvature of initial configuration
print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name);

% Prints .muscle file! [ a_f * Fmax *exp( -( (Q-1)/SK )^2 ) * (1/P0)*(b*P0-a*v)/(v+b); Q = LF/LFO ]
LFO = m_dist; SK = 0.3; a = 0.25; b = 4.0; Fmax = 1e2;
print_Lagrangian_Muscles(xLag,LFO,SK,a,b,Fmax,struct_name,b_ind,t_ind)

% Prints .target file!
%k_Target = 1e7;
%print_Lagrangian_Target_Pts(xLag,k_Target,struct_name);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints VERTEX points to a file called 'struct_name'.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Vertices(xLag,yLag,struct_name)

    N = length(xLag); % Total # of Lag. Pts

    vertex_fid = fopen([struct_name '.vertex'], 'w');

    fprintf(vertex_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        X_v = xLag(s);
        Y_v = yLag(s);
        fprintf(vertex_fid, '%1.16e %1.16e\n', X_v, Y_v);
    end

    fclose(vertex_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints MUSCLE points to a file called "struct_name".muscle
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Muscles(xLag,LFO,SK,a,b,Fmax,struct_name,ind_b,ind_t)

    %N = length(xLag); %Number of Lagrangian Pts. Total

    muscle_fid = fopen([struct_name '.muscle'], 'w');

    fprintf(muscle_fid, '%d\n', 1 ); % Only 1 muscle

    %MUSCLES BETWEEN VERTICES
    fprintf(muscle_fid, '%d %d %1.16e %1.16e %1.16e %1.16e %1.16e\n', ind_b, ind_t, LFO, SK, a,b,Fmax);

    fclose(muscle_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints SPRING points to a file called 'struct_name'.spring
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Springs(xLag,yLag,k_Spring,ds_Rest,struct_name)

    N = length(xLag);

    spring_fid = fopen([struct_name '.spring'], 'w');

    fprintf(spring_fid, '%d\n', N-1 );    % Print # of springs

    %spring_force = kappa_spring*ds/(ds^2);

    %SPRINGS BETWEEN VERTICES
    for s = 1:N
            if s < N
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds_Rest);
            %else
            %    %Case s=N
            %    fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, 1,   k_Spring, ds_Rest);
            end
    end
    fclose(spring_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints TARGET points to a file called 'struct_name'.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Target_Pts(xLag,k_Target,struct_name)

    N = length(xLag);

    target_fid = fopen([struct_name '.target'], 'w');

    fprintf(target_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    fclose(target_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints BEAM (Torsional Spring) points to a file called 'struct_name'.beam
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name)

    % k_Beam: beam stiffness
    % C: beam curvature

    N = length(xLag); % NOTE: Total number of beams = Number of Total Lag Pts. - 2

    beam_fid = fopen([struct_name '.beam'], 'w');

    fprintf(beam_fid, '%d\n', N-2 );

    %spring_force = kappa_spring*ds/(ds^2);

    %BEAMS BETWEEN VERTICES
    for s = 2:N-1
            if  s <= N-1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1, k_Beam, C(s) );
            else
                %Case s=N
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, 1,   k_Beam, C(s) );
            end
    end
    fclose(beam_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: computes "curvature" of starting configuration
%
% NOTE: not curvature in the traditional geometric sense, in the 'discrete'
% sense through cross product.
%
% NOTE: assumes a CLOSED structure
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function C = compute_Curvatures(xLag,yLag)

N = length(xLag);
C = zeros( N );

%Note: needs to be done same order as you print .beam file!
for i=1:N

   % Pts Xp -> Xq -> Xr (same as beam force calc.)

   if ( (i > 1) && (i < N) )

        Xp = xLag(i-1); Xq = xLag(i); Xr = xLag(i+1);
        Yp = yLag(i-1); Yq = yLag(i); Yr = yLag(i+1);

   elseif (i==1)

        Xp = xLag(N); Xq = xLag(i); Xr = xLag(i+1);
        Yp = yLag(N); Yq = yLag(i); Yr = yLag(i+1);

   elseif (i==N)

        Xp = xLag(N-1); Xq = xLag(N); Xr = xLag(1);
        Yp = yLag(N-1); Yq = yLag(N); Yr = yLag(1);

   end

   C(i) = (Xr-Xq)*(Yq-Yp) - (Yr-Yq)*(Xq-Xp); %Cross product btwn vectors

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(ds,Lx)

% ds: Lagrangian pt. spacing
% Lx: Length of Eulerian grid

% The immsersed structure is initially "an angle" of 10 degrees %
ang = 70*pi/180; % original angle
L = 1/10*Lx;     % Structure will be 1/20 of the length of the grid

% Create line of points
x = 0:ds:L;
y = zeros( 1, length( x ) );

% Rotate line of points
[xB,yB] = rotate_geometry_please(x,y,-ang/2);
[xT,yT] = rotate_geometry_please(x,y,ang/2);

% Put together!
xLag = [ xB(end:-1:2) xT ];
yLag = [ yB(end:-1:2) yT ];

[xLag,yLag] = rotate_geometry_please(xLag,yLag,-pi/4);

length(xLag)

xLag = xLag + 0.75*Lx;
yLag = yLag + 0.25*Lx;

% TESTING GEOMETRY
plot(xB,yB,'r*'); hold on;
plot(xT,yT,'go'); hold on;
plot( xLag, yLag, 'k+'); hold on;
axis([0 Lx 0 Lx]);

% for i=1:length(xLag)
%    plot( xLag(i), yLag(i), 'k+'); hold on;
%    axis([0 Lx 0 Lx]);
%    pause();
% end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: rotate set of points
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xN,yN] = rotate_geometry_please(x,y,ang)

% x,y original (x,y) coordinates before rotation
% ang: angle rotating by

xN = zeros( 1, length( x ) ); yN = xN;
for i=1:length(x)
   xN(i) = x(i)*cos(ang) - y(i)*sin(ang);
   yN(i) = x(i)*sin(ang) + y(i)*cos(ang);
end

author: nickabattista