text:
function Stat=clusterstat(S,partition,between)
%function Stat=clusterstat(S,partition,[between])
%
%PURPOSE
%
%To compute various intra- and extra-cluster statistics.
%
%INPUT
%
%[An argument in brackets is optional. If it isn't  given or it's
% an empty matrix/string, the function will use a default value.]
%
% S         (matrix) NxN similarity or distance matrix between N
%             objects. Must be symmetric.
% partition (1xN vector) partition vector (see explanation for
%             'partition vector' in function hcluster)
% [between] (scalar) 1: compute between-cluster similarity/distance matrix,
%                    0: (default) don't compute
%
%OUTPUT
%
% All fields are vectors of size 1xK, where K is the number of clusters
%
% Stat.N(i)             the number of objects in cluster i
%
% Stat.internal.min(i)  minimum/average/max internal similarity/distance
% Stat.internal.avg(i)  in cluster i
% Stat.internal.max(i)
%  Note: if there is only one item in the cluster, the internal statistics
%  are set to NaN
%
% Stat.external.min(i)  minimum/average/max external similarity/distance from
% Stat.external.avg(i)  objects of cluster i to the objects of other clusters
% Stat.external.max(i)
%
%DETAILS
%
%Partition vector divides items  into clusters: partition(i) is
%the label  of cluster that item i belongs to. Cluster labels must
%be integers 1,2,...,K where K is the number of clusters.
%
%"Internal" for cluster k refers to all pairwise
%distances/similarities between objects belonging to the same
%cluster, i.e., for objects i for which partition(i)==k.
%
%"External" for cluster k refers to all pairwise
%distances/similarities between the members of the cluster k and
%members of the other clusters, i.e., the whole similarity matrix
%excluding the internal similarities of cluster k.
%
%If 'between' is set to 1 the function computes also the following
%KxK matrices
%
%Stat.between.min(i,j)
%Stat.between.avg(i,j)
%Stat.between.max(i,j)
%
%These include distances/similarities  between clusters i and j
%defined as min/average/max pairwise distances between objects
%belonging to clusters i and j, respectively.
%
%SEE ALSO
% clusterquality
% icassoStability
% icassoRindex
% rindex

%COPYRIGHT NOTICE
%This function is a part of Icasso software library
%Copyright (C) 2003-2005 Johan Himberg
%
%This program is free software; you can redistribute it and/or
%modify it under the terms of the GNU General Public License
%as published by the Free Software Foundation; either version 2
%of the License, or any later version.
%
%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.
%
%You should have received a copy of the GNU General Public License
%along with this program; if not, write to the Free Software
%Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

% ver 1.2. 100105

if nargin<2,
   error('You must give at least two input arguments.');
end

if nargin<3|isempty(between)
  between=0;
end

% Number of clusters
Ncluster=max(partition);

%Initialize the struct
Stat.internal.sum(1:Ncluster,1)=NaN;
Stat.internal.min(1:Ncluster,1)=NaN;
Stat.internal.avg(1:Ncluster,1)=NaN;
Stat.internal.max(1:Ncluster,1)=NaN;
Stat.external.sum(1:Ncluster,1)=NaN;
Stat.external.min(1:Ncluster,1)=NaN;
Stat.external.avg(1:Ncluster,1)=NaN;
Stat.external.max(1:Ncluster,1)=NaN;
Stat.index = cell(1,Ncluster);

for cluster=1:Ncluster,
  thisPartition=(partition==cluster);
  Stat.index{cluster} = find(thisPartition);
  S_=S(thisPartition,thisPartition);
  Stat.N(cluster)=size(S_,1);
  S_(eye(size(S_))==1)=[];
  if ~isempty(S_),
    Stat.internal.sum(cluster)=sum(S_);
    Stat.internal.min(cluster)=min(S_);
    Stat.internal.avg(cluster)=mean(S_);
    Stat.internal.max(cluster)=max(S_);
  end
  if Ncluster>1,
    S_=S(thisPartition,~thisPartition);
    Stat.external.sum(cluster)=sum(S_(:));
    Stat.external.min(cluster)=min(S_(:));
    Stat.external.avg(cluster)=mean(S_(:));
    Stat.external.max(cluster)=max(S_(:));
  end
end

if between,
  Stat.between.min=zeros(Ncluster,Ncluster);
  Stat.between.max=zeros(Ncluster,Ncluster);
  Stat.between.avg=zeros(Ncluster,Ncluster);

  for i=1:Ncluster,
    Pi=find(i==partition);
    for j=i+1:Ncluster,
      Pj=find(j==partition);
      d_=S(Pi,Pj);
      Stat.between.min(i,j)=min(d_(:));
      Stat.between.avg(i,j)=mean(d_(:));
      Stat.between.max(i,j)=max(d_(:));		;
    end
  end

  Stat.between.min=Stat.between.min+ ...
    Stat.between.min';
  Stat.between.max=Stat.between.max+ ...
      Stat.between.max';
  Stat.between.avg=Stat.between.avg+ ...
      Stat.between.avg';
end

author: fieldtrip

repo: fieldtrip


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: May 27th, 2015
% 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 JELLYFISH-EXAMPLE geometry and prints associated input files
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Jellyfish_Geometry()

% FLUID GRID PARAMETERS %
L = 8;                              % Length of computational domain (m)
N = 256;                            % # of Cartesian grid meshwidths
dx = L/N;                           % Cartesian mesh width (m)
ds = L/(2.0*N);                     % Ideal Lagrangian spacing (m)

% Construct Geometry
[xLag,yLag,ds] = give_Me_Immsersed_Boundary_Geometry(N,L);
plot(xLag,yLag,'*'); hold on;

% Translate Geometry
xLag = xLag + L/8;
yLag = yLag + L/4;
plot(xLag,yLag,'r*'); hold on;

% NAMING CONVENTION FOR SIMULATION
struct_name = 'jelly';      % structure name

%
% PRINT INPUT FILES (.vertex, .spring, .beam, etc) %
%

% print vertices
print_Lagrangian_Vertices(xLag,yLag,struct_name);

% print springs
k_Spring = 500*1.2750000000000000e+07; %500   % spring constant (Newton) dt=2.5e-6
print_Lagrangian_Springs(xLag,yLag,k_Spring,ds,struct_name);

% print beams
k_Beam = 1.0363359375000002e+13;  %100 % beam stiffness constant (Newton m^2) %5.1816796875000010e+12
print_Lagrangian_nonInv_Beams(xLag,yLag,k_Beam,struct_name);

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

function print_Lagrangian_Vertices(xLag,yLag,struct_name)

    N = length(xLag);

    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 SPRING points to a file called rubberband.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 );

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

    % SPRINGS BETWEEN VERTICES ON RHS
    for s = 1:ceil(N/2)
            if s <= floor(N/2)
                x1 = xLag(s); x2 = xLag(s+1);
                y1 = yLag(s); y2 = yLag(s+1);
                ds = sqrt( (x1-x2)^2 + (y1-y2)^2 );
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds);
            elseif s == ceil(N/2)
                x1 = xLag(1); x2 = xLag( ceil(N/2)+1 );
                y1 = yLag(1); y2 = yLag( ceil(N/2)+1 );
                ds = sqrt( (x1-x2)^2 + (y1-y2)^2 );
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', 1, ceil(N/2)+1, k_Spring, ds);
            end
    end

    % SPRINGS BETWEEN VERTICES ON LHS
    for s=1:floor(N/2)-1
                s1 = ceil(N/2)+s;
                s2 = ceil(N/2)+s+1;
                x1 = xLag(s1); x2 = xLag(s2);
                y1 = yLag(s1); y2 = yLag(s2);
                ds = sqrt( (x1-x2)^2 + (y1-y2)^2 );
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s1, s2, k_Spring, ds);
    end
    fclose(spring_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints BEAM (NON-INVARIANT) points to a file called rubberband.nonInv_beam
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_nonInv_Beams(xLag,yLag,k_Beam,struct_name)

    % k_Beam: beam stiffness
    % Cx/Cy:  beam curvatures in x/y respestively

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

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

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

    % beam_force = kappa_beam*ds/(ds^4)

    %BEAMS BETWEEN VERTICES ON RHS
    for s = 1:ceil(N/2)+1
            if s <= floor(N/2)-1
                Cx = xLag(s) - 2*xLag(s+1) + xLag(s+2);
                Cy = yLag(s) - 2*yLag(s+1) + yLag(s+2);
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e %1.16e\n', s,s+1,s+2, k_Beam, Cx, Cy);
            elseif s == ceil(N/2)
                Cx = xLag(ceil(N/2)+1) - 2*xLag(1) + xLag(2);
                Cy = yLag(ceil(N/2)+1) - 2*yLag(1) + yLag(2);
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e %1.16e\n', ceil(N/2)+1,1,2, k_Beam, Cx,Cy);
            elseif s == ceil(N/2)+1
                Cx = xLag(ceil(N/2)+2) - 2*xLag( ceil(N/2)+1 ) + xLag(1);
                Cy = yLag(ceil(N/2)+2) - 2*yLag( ceil(N/2)+1 ) + yLag(1);
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e %1.16e\n', ceil(N/2)+2,ceil(N/2)+1,1, k_Beam, Cx,Cy);
            end
    end

    % BEAMS BETWEEN VERTICES ON LHS
    for s=1:floor(N/2)-2
                s1 = ceil(N/2)+s;
                s2 = ceil(N/2)+s+1;
                s3 = ceil(N/2)+s+2;
                Cx = xLag(s3) - 2*xLag(s2) + xLag(s1);
                Cy = yLag(s3) - 2*yLag(s2) + yLag(s1);
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e %1.16e\n', s3, s2, s1, k_Beam, Cx, Cy);
    end
    fclose(beam_fid);

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

function [xLag,yLag,ds] = give_Me_Immsersed_Boundary_Geometry(N,L)

% JELLYFISH GEOMETRY PARAMETERS %

bell_length = 2;                    % bell length (m)
npts_bell = ceil(2.0*(bell_length/L)*N); % number of points along the length of the entire bell
npts_circ = 1;                    %number of points along the circumference (=1 for 2D)
npts = npts_bell*npts_circ;	      % total number points
ds = bell_length/(npts_bell-1);   % mesh spacing along the length of bell (m)
Zs = 0;                          %distance to top of bell (m)
xb = zeros(npts);                 %holds resting positions to calculate curvatures
zb=zeros(npts);                   %holds resting positions to calculate curvatures

% Values to describe bell shape given in Alben, Peng and Miller %
betao = 0.5;
betam = 0.3;
to = 0.5;
t=0;
gamma = 1;

% Starting values
z = Zs;
r = 0;
x = 0;

% center line of jelly
xLag(1) = x;
yLag(1) = z;
zl=z;
rl=r;

%right side of bell
for s = 1:(ceil(npts_bell/2)-1)
    beta = betao+(betam-betao)*(t/to)^gamma;
    theta = -1.55*(1-exp(-(s)*ds/beta));
    z = zl + ds*sin(theta);
    r = rl + ds*cos(theta);
    x = r;
    zl=z;
    rl=r;
    xLag(s+1)=x;
    yLag(s+1)=z;
   %fprintf(vertex_fid, '%1.16e %1.16e\n', x, z);
end

% reinitialize values for centerline
z = Zs;
r = 0;
zl=z;
rl=r;

%left side of bell
for s = (ceil(npts_bell/2)):npts_bell-2
    s2=s-(ceil(npts_bell/2)-1);
    beta = betao+(betam-betao)*(t/to)^gamma;
    theta = -1.55*(1-exp(-(s2)*ds/beta));
    z = zl + ds*sin(theta);
    r = rl + ds*cos(theta);
    x = -1*r;
    zl=z;
    rl=r;
    xLag(s+1)=x;
    yLag(s+1)=z;
   %fprintf(vertex_fid, '%1.16e %1.16e\n', x, z);
end

%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% % Write out the beam information
%  beam_fid = fopen([mesh_name num2str(N) '.beam'], 'w');
%
%  fprintf(beam_fid, '%d\n', npts-3);
%
%  %right side of bell
%  for q = 0:npts_circ-1
%  for s = 0:(ceil(npts_bell/2)-3)
%     C1 = xb(s+1)+xb(s+3)-2*xb(s+2);
%     C2 = zb(s+1)+zb(s+3)-2*zb(s+2);
%     fprintf(beam_fid, '%d %d %d %1.16e %1.16e %1.16e\n', q*npts_circ+s, q*npts_circ+s+1, q*npts_circ+s+2, kappa_beam*ds/(ds^4), C1, C2);
%  end
%
%  %top of bell
%     s=ceil(npts_bell/2);
%     C1 = xb(s+1)+xb(2)-2*xb(1);
%     C2 = zb(s+1)+zb(2)-2*zb(1);
%     fprintf(beam_fid, '%d %d %d %1.16e %1.16e %1.16e\n', q*npts_circ+s, q*npts_circ+0, q*npts_circ+1, kappa_beam*ds/(ds^4), C1, C2);
%     s=ceil(npts_bell/2)+1;
%     C1 = xb(s+1)+xb(1)-2*xb(s);
%     C2 = zb(s+1)+zb(1)-2*zb(s);
%     fprintf(beam_fid, '%d %d %d %1.16e %1.16e %1.16e\n', q*npts_circ+s, q*npts_circ+s-1, q*npts_circ+0, kappa_beam*ds/(ds^4), C1, C2);
%
%  %left side of bell
%  for s = ceil((npts_bell/2)+2):npts_bell-2
%     C1 = xb(s+1)+xb(s-1)-2*xb(s);
%     C2 = zb(s+1)+zb(s-1)-2*zb(s);
%     fprintf(beam_fid, '%d %d %d %1.16e %1.16e %1.16e\n', q*npts_circ+s, q*npts_circ+s-1, q*npts_circ+s-2, kappa_beam*ds/(ds^4), C1, C2);
%  end
%  end
% %
% %
%  fclose(beam_fid);
%

author: nickabattista


text:
function inside = polygon_contains_point_2d ( n, v, p )

%*****************************************************************************80
%
%% POLYGON_CONTAINS_POINT_2D finds if a point is inside a simple polygon in 2D.
%
%  Discussion:
%
%    A simple polygon is one whose boundary never crosses itself.
%    The polygon does not need to be convex.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    M Shimrat,
%    Position of Point Relative to Polygon,
%    ACM Algorithm 112,
%    Communications of the ACM,
%    Volume 5, Number 8, page 434, August 1962.
%
%  Parameters:
%
%    Input, integer N, the number of nodes or vertices in the polygon.
%    N must be at least 3.
%
%    Input, real V(2,N), the coordinates of the vertices of the polygon.
%
%    Input, real P(2,1), the coordinates of the point to be tested.
%
%    Output, logical INSIDE, is TRUE if the point is inside the polygon.
%
  inside = 0;

  for i = 1 : n

    x1 = v(1,i);
    y1 = v(2,i);

    if ( i < n )
      x2 = v(1,i+1);
      y2 = v(2,i+1);
    else
      x2 = v(1,1);
      y2 = v(2,1);
    end

    if ( ( y1 < p(2,1) & p(2,1) <= y2 ) | ( p(2,1) <= y1 & y2 < p(2,1) ) )
      if ( ( p(1,1) - x1 ) - ( p(2,1) - y1 ) * ( x2 - x1 ) / ( y2 - y1 ) < 0.0 )
        inside = ~inside;
      end
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function ANew=ellipsCoveringEllipsAndPoint(z,A,p,gammaVal,interiorOpt)
%%ELLIPSCOVERINGELLIPSEANDPOINT Given an ellipse (or ellipsoid) and a
%           point, find the minimum volume ellipse/ellipsoid that is
%           centered at the center of the original ellipse and that covers
%           the original ellipse and the point.
%
%INPUTS: z  The numDimX1 center of the ellipsoid.
%        A  A numDimXnumDim symmetric, positive definite matrix that
%           specifies the size and shape of the ellipse or ellipsoid, where
%           a point zp is on the ellipse/ellipsoid if
%           (zp-z)'*A*(zp-z)=gammaVal.
%        p  A numDimX1 point.
%  gammaVal The threshold for declaring a point to be in the ellipsoid. If
%           this parameter is omitted or an empty matrix is passed, the
%           default value of 1 is used.
% interiorOpt An optional parameter specifying what should be done if the
%           point p is inside of the ellipsoid. Possible values are
%           0 (The default if omitted or an empty matrix is passed) The
%             original ellipsoid is returned if the point is inside it.
%           1 The ellipse that is returned is the maximum volume ellipsoid
%             such that p is on the boundary and the ellipse is entirely
%             contained in the original ellipsoid.
%
%OUTPUTS: ANew The numDimXnumDim new positive definite matrix specifying
%              the shape of the minimum volume ellipsoid covering the
%              original ellipsoid and the point.
%
%The algorithm is taken from Section 10 of [1]. The basic idea is that
%first, we divide A by gammaVal so that the equation for the ellipse is
%(zp-z)'*A*(zp-z)=1. Next, let L be the lower-triangular Cholesky
%decomposition of A. We then do the transformation y=L'*(zp-z) and in the
%transformed domain, the ellipse is the unit sphere norm(y)=1. We can then
%apply the transformation to p to get pTilde=L'*(p-z). We basically now
%just want to transform p so that it is on the x-axis, and then expand the
%sphere in that direction minimally to intersect p and then transform back.
%
%The sphere in y has an equvalent matrix ATilde=eye(numDiom,numDim). We
%want to expand it just enough in a single direction to touch pTilde. This
%is a rank-one modification problem, which is described in a little detail
%in Sections 4 and 10 of [1].
%
%EXAMPLE 1:
%Here we have an ellipse and a point outside of the ellipse.
% A=[1,0;
%    0,10];
% M=[0.413074198133900,  0.910697373904216;
%    -0.910697373904216,   0.413074198133900];%A rotation matrix
% A=M*A*M';
% z=[5;6];
% p=[6;7];
% ANew=ellipsCoveringEllipsAndPoint(z,A,p);
% figure(1)
% clf
% hold on
% drawEllipse(z,A,1,'b','linewidth',2)
% drawEllipse(z,ANew,1,'m','linewidth',2)
% scatter(p(1),p(2),'xr','linewidth',2)
%The new ellipse engulf the original ellipse and the point.
%
%EXAMPLE 2:
%In this instance, we have an ellipse and a point inside the ellipse.
% A=[1,0;
%    0,10];
% M=[0.413074198133900,  0.910697373904216;
%    -0.910697373904216,   0.413074198133900];%A rotation matrix
% A=M*A*M';
% z=[5;6];
% p=[5.2;6];
% ANew=ellipsCoveringEllipsAndPoint(z,A,p);
% figure(1)
% clf
% hold on
% drawEllipse(z,A,1,'b','linewidth',2)
% drawEllipse(z,ANew,1,'m','linewidth',2)
% scatter(p(1),p(2),'xr','linewidth',2)
%The new ellipse is the same as the orignal ellipse, because the point is
%inside the original ellipse.
%
%REFERENCES:
%[1] S. B. Pope, "Algorithms for ellipsoids," Cornell University, Tech.
%    Rep. FDA-08-01, Feb. 2008. [Online].
%    Available: https://tcg.mae.cornell.edu/pubs/Pope_FDA_08.pdf
%
%May 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<5||isempty(interiorOpt))
    interiorOpt=0;
end

if(nargin<4||isempty(gammaVal))
   gammaVal=1;
end

n=length(z);

A=A/gammaVal;

%Check for the point being inside the ellipse.
diff=p-z;
if(interiorOpt==0&&diff'*A*diff<=1)
    %The point is already inside the ellipse.
    ANew=A;
    return;
end

L=chol(A,'lower');

pTilde=L'*(p-z);%Eq 62.

pN=norm(pTilde);
g=(1/pN-1)*(1/pN^2);
G=eye(n,n)+g*(pTilde*pTilde');
ANew=L*G*(L*G)'*gammaVal;

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function u = sqr(a)
%SQR          Slope square  sqr(a)
%

% written  12/06/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  INTLAB_SLOPE = getappdata(0,'INTLAB_SLOPE');

  u = a;

  u.r = sqr(a.r);
  indexc = 1:INTLAB_SLOPE.NUMVAR;
  indexr = 2:INTLAB_SLOPE.NUMVAR+1;
  u.s = a.s .* (a.r(:,indexc)+a.r(:,indexr));

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function [ Y ] = AfuncAcc( U1, V1, U0, V0, spa, R)

Y = U1*(V1'*R);
Y = Y + U0*(V0'*R);
Y = Y + spa*R;

end

author: HKUST-KnowComp

repo: FMG

save_path: github-repos/MATLAB/HKUST-KnowComp-FMG


text:
function exact = p35_exact ( )

%*****************************************************************************80
%
%% P35_EXACT returns the exact integral for problem 35.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the value of the integral.
%
  exact = 26.0;

  return
end


text:
% tidy [(position)]
% arranges the command window, editor window, and other specified windows
% neatly on screen.  easily customised by editing script.
% choose the command window position from
%   tidy bottom
%   tidy top
%   tidy left
%   tidy right
%   tidy compact
%   tidy (defaults to bottom)
% should work under win95/98/NT/2000, tested on win2000 only
% should work using any version of MATLAB, tested using MATLAB5/6 only

function tidy(pos,extra)

% also can be used as a shortcut to MoveWindow for special cases
if nargin==2
    MoveWindow(pos,extra)
    return;
end

if nargin==0, pos='bottom'; end

switch pos
case {'top','bottom','left','right'}
    disp('tidying desktop...');
    compact=0;
case {'compact'}
    disp('tidying desktop...');
    pos='bottom';
    compact=1;
otherwise
    help tidy
    disp(['Unrecognised position "' pos '".']);
    return;
end

% all windows starting with the supplied title are matched, but
% typically only the first one encountered is moved.
% a '*' as the last character of a window title means move all
% of the windows that match.  typically this might be used
% as below - 'Figure No. *' will move all figure windows that
% have not had their captions altered.
% the command window should probably be last in the list since
% whichever is last in the list has the focus after these
% operations.
MoveWindow('Figure No. *',SecondDivision(Opposite(pos)));
MoveWindow('Figure No. 2',FirstDivision(Opposite(pos)));
MoveWindow('Microsoft Visual C++',Opposite(pos));
MoveWindow('D:\*',SecondDivision(Opposite(pos)));
MoveWindow('d:\*',SecondDivision(Opposite(pos)));
MoveWindow('E:\*',FirstDivision(Opposite(pos)));
MoveWindow('e:\*',FirstDivision(Opposite(pos)));
MoveWindow('Workspace',SecondDivision(Opposite(pos)));
if compact
    MoveWindow('MATLAB',FirstDivision(pos));
    MoveWindow('ALTOUT',SecondDivision(pos));
else
    MoveWindow('MATLAB',pos);
    MoveWindow('ALTOUT',SecondDivision(Opposite(pos)));
end

function opp=Opposite(pos)
if strcmp(pos,'top') opp='bottom'; end
if strcmp(pos,'bottom') opp='top'; end
if strcmp(pos,'left') opp='right'; end
if strcmp(pos,'right') opp='left'; end
if strcmp(pos,'topleft') opp='bottomright'; end
if strcmp(pos,'topright') opp='bottomleft'; end
if strcmp(pos,'bottomleft') opp='topright'; end
if strcmp(pos,'bottomright') opp='topleft'; end

function div=FirstDivision(pos)
if strcmp(pos,'top') div='topleft'; end
if strcmp(pos,'bottom') div='bottomleft'; end
if strcmp(pos,'left') div='topleft'; end
if strcmp(pos,'right') div='topright'; end

function div=SecondDivision(pos)
if strcmp(pos,'top') div='topright'; end
if strcmp(pos,'bottom') div='bottomright'; end
if strcmp(pos,'left') div='bottomleft'; end
if strcmp(pos,'right') div='bottomright'; end

author: Sable

repo: mcbench-benchmarks


text:
function [Q_matrix] = adapt_noise_covariance(Phi_matrix, P_new, P_old, Q_matrix, n, GNSS_epoch, corrections)
%% Adapt Q matrix to system noise
% Adjust measurement noise covariance matrix Q according to state correction sequence
% Adam Werries 2016, see Apache 2.0 license.

%% Mohamed and Schwarz method, 1999
k = GNSS_epoch-1;
if k > n
    C = zeros(15,15);
    for j = k-n+1:k
        C = C + corrections(:,j)*corrections(:,j)';
    end
    Q_matrix = C./n + P_new - Phi_matrix*P_old*Phi_matrix';
    Q_matrix = diag(diag(Q_matrix));
else

end

%% Ding and Wang method, 2007
% TODO

end

author: awerries

repo: kalman-localization


text:
filename='Bridge_hexahedra_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'SLERP'; incrementFactor = 1;
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.05;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

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

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [V,nr,nre]=lcon2vert(A,b,Aeq,beq,TOL,checkbounds)
import iris.thirdParty.polytopes.*;
%An extension of Michael Kleder's con2vert function, used for finding the
%vertices of a bounded polyhedron in R^n, given its representation as a set
%of linear constraints. This wrapper extends the capabilities of con2vert to
%also handle cases where the  polyhedron is not solid in R^n, i.e., where the
%polyhedron is defined by both equality and inequality constraints.
%
%SYNTAX:
%
%  [V,nr,nre]=lcon2vert(A,b,Aeq,beq,TOL)
%
%The rows of the N x n matrix V are a series of N vertices of the polyhedron
%in R^n, defined by the linear constraints
%
%   A*x  <= b
%   Aeq*x = beq
%
%By default, Aeq=beq=[], implying no equality constraints. The output "nr"
%lists non-redundant inequality constraints, and "nre" lists non-redundant
%equality constraints.
%
%The optional TOL argument is a tolerance used for both rank-estimation and
%for testing feasibility of the equality constraints. Default=1e-10.
%The default can also be obtained by passing TOL=[];
%
%
%EXAMPLE:
%
%The 3D region defined by x+y+z=1, x>=0, y>=0, z>=0
%is described by the following constraint data.
%
%
%     A =
%
%         0.4082   -0.8165    0.4082
%         0.4082    0.4082   -0.8165
%        -0.8165    0.4082    0.4082
%
%
%     b =
%
%         0.4082
%         0.4082
%         0.4082
%
%
%     Aeq =
%
%         0.5774    0.5774    0.5774
%
%
%     beq =
%
%         0.5774
%
%
%  >> V=lcon2vert(A,b,Aeq,beq)
%
%         V =
%
%             1.0000    0.0000    0.0000
%             0.0000    0.0000    1.0000
%            -0.0000    1.0000    0.0000
%
%

  %%initial argument parsing

  nre=[];
  nr=[];
  if nargin<5 || isempty(TOL), TOL=1e-10; end
  if nargin<6, checkbounds=true; end

  switch nargin

      case 0

           error 'At least 1 input argument required'

      case 1

         b=[]; Aeq=[]; beq=[];

      case 2

          Aeq=[]; beq=[];

      case 3

          beq=[];
          error 'Since argument Aeq specified, beq must also be specified'

  end

  b=b(:); beq=beq(:);

  if xor(isempty(A), isempty(b))
     error 'Since argument A specified, b must also be specified'
  end

  if xor(isempty(Aeq), isempty(beq))
        error 'Since argument Aeq specified, beq must also be specified'
  end

  nn=max(size(A,2)*~isempty(A),size(Aeq,2)*~isempty(Aeq));

  if ~isempty(A) && ~isempty(Aeq) && ( size(A,2)~=nn || size(Aeq,2)~=nn)

      error 'A and Aeq must have the same number of columns if both non-empty'

  end

  inequalityConstrained=~isempty(A);
  equalityConstrained=~isempty(Aeq);

 [A,b]=rownormalize(A,b);
 [Aeq,beq]=rownormalize(Aeq,beq);

  if equalityConstrained && nargout>2

        [discard,nre]=lindep([Aeq,beq].',TOL);

        if ~isempty(nre) %reduce the equality constraints

            Aeq=Aeq(nre,:);
            beq=beq(nre);

        else
            equalityConstrained=false;
        end

   end

   %%Find 1 solution to equality constraints within tolerance

   if equalityConstrained

       Neq=null(Aeq);

       x0=pinv(Aeq)*beq;

       if norm(Aeq*x0-beq)>TOL*norm(beq),  %infeasible

          nre=[]; nr=[]; %All constraints redundant for empty polytopes
          V=[];
          return;

       elseif isempty(Neq)

           V=x0(:).';
           nre=(1:nn).'; %Equality constraints determine everything.
           nr=[];%All inequality constraints are therefore redundant.
           return

       end

       rkAeq= nn - size(Neq,2);

  end

    %%
  if inequalityConstrained && equalityConstrained

   AAA=A*Neq;
   bbb=b-A*x0;

  elseif inequalityConstrained

    AAA=A;
    bbb=b;

  elseif equalityConstrained && ~inequalityConstrained

       error('Non-bounding constraints detected. (Consider box constraints on variables.)')

  end

  nnn=size(AAA,2);

  if nnn==1 %Special case

     idxu=sign(AAA)==1;
     idxl=sign(AAA)==-1;
     idx0=sign(AAA)==0;

     Q=bbb./AAA;
     U=Q;
       U(~idxu)=inf;
     L=Q;
       L(~idxl)=-inf;

     [ub,uloc]=min(U);
     [lb,lloc]=max(L);

     if ~all(bbb(idx0)>=0) || ub<lb %infeasible

         V=[]; nr=[]; nre=[];
         return

     elseif ~isfinite(ub) || ~isfinite(lb)

         error('Non-bounding constraints detected. (Consider box constraints on variables.)')

     end

     Zt=[lb;ub];

     if nargout>1
        nr=unique([lloc,uloc]); nr=nr(:);
     end

  else

          if nargout>1
           [Zt,nr]=con2vert(AAA,bbb,TOL,checkbounds);
          else
            Zt=con2vert(AAA,bbb,TOL,checkbounds);
          end

  end

  if equalityConstrained && ~isempty(Zt)

      V=bsxfun(@plus,Zt*Neq.',x0(:).');

  else

      V=Zt;

  end

  if isempty(V),
     nr=[]; nre=[];
  end

 function [V,nr] = con2vert(A,b,TOL,checkbounds)
% CON2VERT - convert a convex set of constraint inequalities into the set
%            of vertices at the intersections of those inequalities;i.e.,
%            solve the "vertex enumeration" problem. Additionally,
%            identify redundant entries in the list of inequalities.
%
% V = con2vert(A,b)
% [V,nr] = con2vert(A,b)
%
% Converts the polytope (convex polygon, polyhedron, etc.) defined by the
% system of inequalities A*x <= b into a list of vertices V. Each ROW
% of V is a vertex. For n variables:
% A = m x n matrix, where m >= n (m constraints, n variables)
% b = m x 1 vector (m constraints)
% V = p x n matrix (p vertices, n variables)
% nr = list of the rows in A which are NOT redundant constraints
%
% NOTES: (1) This program employs a primal-dual polytope method.
%        (2) In dimensions higher than 2, redundant vertices can
%            appear using this method. This program detects redundancies
%            at up to 6 digits of precision, then returns the
%            unique vertices.
%        (3) Non-bounding constraints give erroneous results; therefore,
%            the program detects non-bounding constraints and returns
%            an error. You may wish to implement large "box" constraints
%            on your variables if you need to induce bounding. For example,
%            if x is a person's height in feet, the box constraint
%            -1 <= x <= 1000 would be a reasonable choice to induce
%            boundedness, since no possible solution for x would be
%            prohibited by the bounding box.
%        (4) This program requires that the feasible region have some
%            finite extent in all dimensions. For example, the feasible
%            region cannot be a line segment in 2-D space, or a plane
%            in 3-D space.
%        (5) At least two dimensions are required.
%        (6) See companion function VERT2CON.
%        (7) ver 1.0: initial version, June 2005
%        (8) ver 1.1: enhanced redundancy checks, July 2005
%        (9) Written by Michael Kleder
%
%Modified by Matt Jacobson - March 30, 2011
%
  import iris.thirdParty.polytopes.*;

   %%%3/4/2012 Improved boundedness test - unfortunately slower than Michael Kleder's
   if checkbounds

    [aa,bb,aaeq,bbeq]=vert2lcon(A,TOL);

    if any(bb<=0) || ~isempty(bbeq)
        error('Non-bounding constraints detected. (Consider box constraints on variables.)')
    end

    clear aa bb aaeq bbeq

   end

   dim=size(A,2);

   %%%Matt J initialization
   if strictinpoly(b,TOL)

       c=zeros(dim,1);

   else

            slackfun=@(c)b-A*c;

            %Initializer0
            c = pinv(A)*b; %02/17/2012 -replaced with pinv()
            s=slackfun(c);

            if ~approxinpoly(s,TOL) %Initializer1

                c=Initializer1(TOL,A,b,c);
                s=slackfun(c);

            end

            if  ~approxinpoly(s,TOL)  %Attempt refinement

                %disp 'It is unusually difficult to find an interior point of your polytope. This may take some time... '
                %disp ' '

                c=Initializer2(TOL,A,b,c);
                %[c,fval]=Initializer1(TOL,A,b,c,10000);
                s=slackfun(c);

            end

            if ~approxinpoly(s,TOL)
                    %error('Unable to locate a point near the interior of the feasible region.')
                    V=[];
                    nr=[];
                    return
            end

           if ~strictinpoly(s,TOL) %Added 02/17/2012 to handle initializers too close to polytope surface

                %disp 'Recursing...'

                idx=(  abs(s)<=max(s)*TOL );

                Amod=A; bmod=b;
                 Amod(idx,:)=[];
                 bmod(idx)=[];

                Aeq=A(idx,:); %pick the nearest face to c
                beq=b(idx);

                faceVertices=lcon2vert(Amod,bmod,Aeq,beq,TOL,1);
                if isempty(faceVertices)
                   disp 'Something''s wrong. Couldn''t find face vertices. Possibly polyhedron is unbounded.'
                   keyboard
                end

                c=faceVertices(1,:).';  %Take any vertex - find local recession cone vector
                s=slackfun(c);

                idx=(  abs(s)<=max(s)*TOL );

                Asub=A(idx,:); bsub=b(idx,:);

                [aa,bb,aaeq,bbeq]=vert2lcon(Asub);
                aa=[aa;aaeq;-aaeq];
                bb=[bb;bbeq;-bbeq];

                clear aaeq bbeq

                [bmin,idx]=min(bb);

                 if bmin>=-TOL
                   disp 'Something''s wrong. We should have found a recession vector (bb<0).'
                   keyboard
                 end

                Aeq2=null(aa(idx,:)).';
                beq2=Aeq2*c;  %find intersection of polytope with line through facet centroid.

                linetips = lcon2vert(A,b,Aeq2,beq2,TOL,1);

                if size(linetips,1)<2
                   disp 'Failed to identify line segment through interior.'
                   disp 'Possibly {x: Aeq*x=beq} has weak intersection with interior({x: Ax<=b}).'
                   keyboard
                end

                lineCentroid=mean(linetips);%Relies on boundedness

                clear aa bb

                c=lineCentroid(:);
                s=slackfun(c);

            end

            b = s;
   end
    %%%end Matt J initialization

    D=bsxfun(@rdivide,A,b);

    k = convhulln(D);
    nr = unique(k(:));

    G  = zeros(size(k,1),dim);
    ee=ones(size(k,2),1);
    discard=false( 1, size(k,1) );

    for ix = 1:size(k,1) %02/17/2012 - modified

        F = D(k(ix,:),:);
        if lindep(F,TOL)<dim;
            discard(ix)=1;
            continue;
        end

        G(ix,:)=F\ee;

    end

    G(discard,:)=[];

    V = bsxfun(@plus, G, c.');

    [discard,I]=unique( round(V*1e6),'rows');
    V=V(I,:);

return

function [c,fval]=Initializer1(TOL, A,b,c,maxIter)

    thresh=-10*max(eps(b));

    if nargin>4
     [c,fval]=fminsearch(@(x) max([thresh;A*x-b]), c,optimset('MaxIter',maxIter));
    else
     [c,fval]=fminsearch(@(x) max([thresh;A*x-b]), c);
    end

return

function c=Initializer2(TOL,A,b,c)
 %norm(  (I-A*pinv(A))*(s-b) )  subj. to s>=0

    maxIter=100000;

    [mm,nn]=size(A);

     Ap=pinv(A);
     Aaug=speye(mm)-A*Ap;
     Aaugt=Aaug.';

    M=Aaugt*Aaug;
    C=sum(abs(M),2);
     C(C<=0)=min(C(C>0));

    slack=b-A*c;
    slack(slack<0)=0;

        %     relto=norm(b);
        %     relto =relto + (relto==0);
        %
        %      relres=norm(A*c-b)/relto;

    IterThresh=maxIter;
    s=slack;
    ii=0;
    %for ii=1:maxIter
    while ii<=2*maxIter %HARDCODE

       ii=ii+1;
       if ii>IterThresh,
           %warning 'This is taking a lot of iterations'
           IterThresh=IterThresh+maxIter;
       end

     s=s-Aaugt*(Aaug*(s-b))./C;
     s(s<0)=0;

       c=Ap*(b-s);
       %slack=b-A*c;
       %relres=norm(slack)/relto;
       %if all(0<slack,1)||relres<1e-6||ii==maxIter, break;  end

    end

return

function [r,idx,Xsub]=lindep(X,tol)
%Extract a linearly independent set of columns of a given matrix X
%
%    [r,idx,Xsub]=lindep(X)
%
%in:
%
%  X: The given input matrix
%  tol: A rank estimation tolerance. Default=1e-10
%
%out:
%
% r: rank estimate
% idx:  Indices (into X) of linearly independent columns
% Xsub: Extracted linearly independent columns of X

   if ~nnz(X) %X has no non-zeros and hence no independent columns

       Xsub=[]; idx=[];
       return
   end

   if nargin<2, tol=1e-10; end

     [Q, R, E] = qr(X,0);

     diagr = abs(diag(R));

     %Rank estimation
     r = find(diagr >= tol*diagr(1), 1, 'last'); %rank estimation

     if nargout>1
      idx=sort(E(1:r));
        idx=idx(:);
     end

     if nargout>2
      Xsub=X(:,idx);
     end

 function [A,b]=rownormalize(A,b)
 %Modifies A,b data pair so that norm of rows of A is either 0 or 1

  if isempty(A), return; end

  normsA=sqrt(sum(A.^2,2));
  idx=normsA>0;
  A(idx,:)=bsxfun(@rdivide,A(idx,:),normsA(idx));
  b(idx)=b(idx)./normsA(idx);

 function tf=approxinpoly(s,TOL)

   smax=max(s);

   if smax<=0
      tf=false; return
   end

   tf=all(s>=-smax*TOL);

  function tf=strictinpoly(s,TOL)

   smax=max(s);

   if smax<=0
      tf=false; return
   end

   tf=all(s>=smax*TOL);

author: rdeits

repo: iris-distro

save_path: github-repos/MATLAB/rdeits-iris-distro


text:
function [ x, y ] = sswap ( n, x, incx, y, incy )

%*****************************************************************************80
%
%% SSWAP interchanges two vectors.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 May 2006
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for Fortran Usage,
%    Algorithm 539,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, September 1979, pages 308-323.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vectors.
%
%    Input, real X(*), one of the vectors to swap.
%
%    Input, integer INCX, the increment between successive entries of X.
%
%    Input, real Y(*), one of the vectors to swap.
%
%    Input, integer INCY, the increment between successive elements of Y.
%
%    Output, real X(*), the swapped vector.
%
%    Output, real Y(*), the swapped vector.
%
  if ( n <= 0 )

  elseif ( incx == 1 & incy == 1 )

    m = mod ( n, 3 );

    for i = 1 : m
      temp = x(i);
      x(i) = y(i);
      y(i) = temp;
    end

    for i = m+1 : 3 : n

      temp = x(i);
      x(i) = y(i);
      y(i) = temp;

      temp = x(i+1);
      x(i+1) = y(i+1);
      y(i+1) = temp;

      temp = x(i+2);
      x(i+2) = y(i+2);
      y(i+2) = temp;

    end

  else

    if ( 0 <= incx )
      ix = 1;
    else
      ix = ( - n + 1 ) * incx + 1;
    end

    if ( 0 <= incy )
      iy = 1;
    else
      iy = ( - n + 1 ) * incy + 1;
    end

    for i = 1 : n
      temp = x(ix);
      x(ix) = y(iy);
      y(iy) = temp;
      ix = ix + incx;
      iy = iy + incy;
    end

  end

  return
end


text:
function p = prior_invunif(varargin)
%PRIOR_INVUNIF  Uniform prior structure for the inverse of the parameter
%
%  Description
%    P = PRIOR_INVUNIF creates uniform prior structure for the
%    inverse of the parameter.
%
%  See also
%    PRIOR_*
%
% Copyright (c) 2009 Jarno Vanhatalo
% Copyright (c) 2010,2012 Aki Vehtari

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

  ip=inputParser;
  ip.FunctionName = 'PRIOR_INVUNIFORM';
  ip.addOptional('p', [], @isstruct);
  ip.parse(varargin{:});
  p=ip.Results.p;

  if isempty(p)
    init=true;
    p.type = 'Inv-Uniform';
  else
    if ~isfield(p,'type') && ~isequal(p.type,'Inv-Uniform')
      error('First argument does not seem to be a valid prior structure')
    end
    init=false;
  end

  if init
    % set functions
    p.fh.pak = @prior_invunif_pak;
    p.fh.unpak = @prior_invunif_unpak;
    p.fh.lp = @prior_invunif_lp;
    p.fh.lpg = @prior_invunif_lpg;
    p.fh.recappend = @prior_invunif_recappend;
  end

end

function [w, s, h] = prior_invunif_pak(p, w)
  w=[];
  s={};
  h=[];
end

function [p, w] = prior_invunif_unpak(p, w)
  w = w;
  p = p;
end

function lp = prior_invunif_lp(x, p)
  lJ=-log(x)*2;   % log(1/x^2) log(|J|) of transformation
  lp = sum(0 +lJ);
end

function lpg = prior_invunif_lpg(x, p)
  lJg=-2./x;      % gradient of log(|J|) of transformation
  lpg = zeros(size(x)) + lJg;
end

function rec = prior_invunif_recappend(rec, ri, p)
% The parameters are not sampled in any case.
  rec = rec;
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function [y, cost] = sig_manifWB(Fopt, FRF, FBB)
[Nt, NRF] = size(FRF);
K = size(FBB,3);

manifold = complexcirclefactory(Nt*NRF);
problem.M = manifold;

for k = 1:K
    temp = Fopt(:,:,k);
    A = kron(FBB(:,:,k).', eye(Nt));
    C1(:,:,k) = temp(:)'*A;
    C2(:,k) = A'*temp(:);
    C3(:,:,k) = A'*A;
    C4(k) = norm(temp,'fro')^2;
end
B1 = sum(C1,3);
B2 = sum(C2,2);
B3 = sum(C3,3);
B4 = sum(C4);

problem.cost = @(x) -B1*x - x'*B2 + trace(B3*x*x') + B4;
problem.egrad = @(x) -2*B2 + 2*B3*x;

% checkgradient(problem);
warning('off', 'manopt:getHessian:approx');

[x,cost] = conjugategradient(problem,FRF(:));
% [x,cost,info,options] = trustregions(problem, FRF(:));
y = reshape(x,Nt,NRF);

end

% problem.cost = @(x) mycost(Fopt, FBB, x);
% function g = mycost(Fopt, FBB, x)
%     g = 0;
%     for k = 1:K
%         temp = Fopt(:,:,k);
%         A = kron(FBB(:,:,k).', eye(Nt));
%         g = g + ( temp(:) -  A*x )' * ( temp(:) - A*x );
%     end
% end
%
% problem.egrad = @(x) mygrad(Fopt, FBB, x);
% function g = mygrad(Fopt, FBB, x)
%     g = 0;
%     for k = 1:K
%         temp = Fopt(:,:,k);
%         A = kron(FBB(:,:,k).', eye(Nt));
%         g = g -2*A'*temp(:) + 2*A'*A*x;
%     end
% end

% problem.costgrad = @(x) mycostgrad(Fopt, FBB, x);
% function [h, g] = mycostgrad(Fopt, FBB, x)
%     h = 0;
%     g = 0;
%     for k = 1:K
%         temp = Fopt(:,:,k);
%         A = kron(FBB(:,:,k).', eye(Nt));
%         g = g -2*A'*temp(:) + 2*A'*A*x;
%         h = h + ( temp(:) -  A*x )' * ( temp(:) - A*x );
%     end
% end

author: Zzhaoxingyu

repo: hybrid-beamforming-for-three-scenes


text:
function fd1d_advection_ftcs ( )

%*****************************************************************************80
%
%% FD1D_ADVECTION_FTCS solves the advection equation using the FTCS method.
%
%  Discussion:
%
%    The FTCS method is unstable for the advection problem.
%
%    Given a smooth initial condition, successive FTCS approximations will
%    exhibit erroneous oscillations of increasing magnitude.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD1D_ADVECTION_FTCS:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solve the constant-velocity advection equation in 1D,\n' );
  fprintf ( 1, '    du/dt = - c du/dx\n' );
  fprintf ( 1, '  over the interval:\n' );
  fprintf ( 1, '    0.0 <= x <= 1.0\n' );
  fprintf ( 1, '  with periodic boundary conditions, and\n' );
  fprintf ( 1, '  with a given initial condition\n' );
  fprintf ( 1, '    u(0,x) = (10x-4)^2 (6-10x)^2 for 0.4 <= x <= 0.6\n' );
  fprintf ( 1, '           = 0 elsewhere.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We use a method known as FTCS:\n' );
  fprintf ( 1, '   FT: Forward Time  : du/dt = (u(t+dt,x)-u(t,x))/dt\n' );
  fprintf ( 1, '   CS: Centered Space: du/dx = (u(t,x+dx)-u(t,x-dx))/2/dx\n' );

  nx = 101;
  dx = 1.0 / ( nx - 1 );
  x = linspace ( 0.0, 1.0, nx );
  nt = 1000;
  dt = 1.0 / nt;
  c = 1.0;

  u = zeros ( 1, nx );
  i = find ( 0.4 <= x & x <= 0.6 );
  u(i) = ( 10.0 * x(i) - 4.0 ).^2 .* ( 6.0 - 10.0 * x(i) ).^2;

  iplot = 1;
  uplot(iplot,:) = u(:)';
  tplot(iplot) = 0.0;
  plotstep = ceil ( nt / 50 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nodes NX = %d\n', nx );
  fprintf ( 1, '  Number of time steps NT = %d\n', nt );
  fprintf ( 1, '  Constant velocity C = %g\n', c );

  im1 = [ nx, 1:nx-2, nx-1 ];
  i   = [ 1,  2:nx-1, nx   ];
  ip1 = [ 2,  3:nx,   1    ];

  for j = 1 : nt

    unew(i) = u(i) - c * dt / dx / 2.0 * ( u(ip1) - u(im1) );
    u(i) = unew(i);

    if ( rem ( j, plotstep ) < 1 )
      iplot = iplot + 1;
      uplot(iplot,:) = u(:)';
      tplot(iplot) = j * dt;
    end

  end
%
%  Plot.
%
  mesh ( x, tplot, uplot );
  xlabel ( '<--X-->' );
  ylabel ( '<--T-->');
  title ( 'U(X,T)');

  filename = 'fd1d_advection_ftcs.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saving plot as "%s".\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD1D_ADVECTION_FTCS\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  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


text:
function varargout = kochsnow(iters,varargin)
% varargout = kochsnow(iters,varargin)
% makes a picture of the Koch snowflake(utilizing draw_bndry and draw_circ)
% made of code from http://www.mathworks.de/matlabcentral/newsreader/view_thread/138798
% credits for the snowflake calculation go to Maarten van Reeuwijk
% varargin examples:
% varargin = 'm', m
% varargin = 'save', 'kochsnow%.bmp'  % is replaced by iters
% varargin = 'show', false
% varargout:
% can be called with one or no argument as an output

    angle = [0, -2/3*pi, -2/3*pi];
    l = 1; m = 1024;
    fileout = false;
    show = true;

    % varargin
    j = 1;
    while j <= nargin-1
        switch varargin{j}
            case 'm'
                m = varargin{j+1}; j = j+2;
            case 'save'
                fileout = varargin{j+1};
                incliters = strfind(fileout,'%');
                sz = size(incliters,2);
                for i = sz:-1:1
                    fileout = strcat(fileout(1:incliters(i)-1), num2str(iters), fileout(incliters(i)+1:end));
                end
                extout = fileout(end-2:end);
                j = j+2;
            case 'show'
                show = varargin{j+1}; j = j+2;
            otherwise
                error('Argument not acceptable!');
        end
    end

    % snowflake calculation by Maarten van Reeuwijk
    for i=1:iters
        l = l/3;
        angle1 = zeros([4*length(angle),1]);
        for j=1:length(angle)
          angle1(4*j-3:4*j) = [angle(j), pi/3, -2*pi/3, pi/3];
        end
        angle = angle1;
    end

    x = zeros([length(angle)+1, 1]);
    y = zeros([length(angle)+1, 1]);

    x(1) = 0; y(1) = 0;
    phi=0;
    for i=1:length(angle);
        phi = phi+angle(i);
        x(i+1) = x(i) + l * cos(phi);
        y(i+1) = y(i) + l * sin(phi);
    end
    % end of snowflake calculation by Maarten van Reeuwijk

    I = false(m*(max(y)-min(y)),m);
    x = (x+min(x))*m; y = (y-min(y))*m;
    I = draw_clline(I, [y,x], 1, 0);
    I = imfill(I,sub2ind(size(I),m/2,m/2),4);

    if show == true, imshow(I,'Border','tight'); end
    if fileout ~= false
        imwrite(I,fileout, extout);
    end
    if nargout == 1, varargout = {I}; end
end


text:
% VL_HOMKERMAP Homogeneous kernel map
%   V = VL_HOMKERMAP(X, N) computes a 2*N+1 dimensional approximated
%   kernel map for the Chi2 kernel. X is an array of data points. Each
%   point is expanded into a vector of dimension 2*N+1 and saved to
%   the output V. The expanded feature vectors are stacked along the
%   first dimension, so that the output array V has the same
%   dimensions of the input array X except for the first one, which is
%   2*N+1 times larger.
%
%   The function accepts the following options:
%
%   Kernel:: KCHI2
%     One of KCHI2 (Chi2 kernel), KINTERS (intersection kernel), KJS
%     (Jensen-Shannon kernel). The 'Kernel' option name can be omitted,
%     i.e. VL_HOMKERMAP(..., 'kernel', 'kchi2') has the same effect of
%     VL_HOMKERMAP(..., 'kchi2').
%
%   Period:: [automatically tuned]
%     Set the period of the kernel specturm. The approximation is
%     based on periodicizing the kernel specturm. If not specified,
%     the period is automatically set based on the heuristic described
%     in [2].
%
%   Window:: [RECTANGULAR]
%     Set the window used to truncate the spectrum before The window
%     can be either RECTANGULAR or UNIFORM window. See [2] and the API
%     documentation for details.
%
%   Gamma:: [1]
%     Set the homogeneity degree of the kernel. The standard kernels
%     are 1-homogeneous, but sometimes smaller values perform better
%     in applications. See [2] for details.
%
%   Example::
%     The following code results in approximatively the same
%     similarities matrices between points X and Y:
%
%       x = rand(10,1) ;
%       y = rand(10,100) ;
%       psix = vl_homkermap(x, 3) ;
%       psiy = vl_homkermap(y, 3) ;
%       figure(1) ; clf ;
%       ker = vl_alldist(x, y, 'kchi2') ;
%       ker_ = psix' * psiy ;
%       plot([ker ; ker_]') ;
%
%   Note::
%     The homogeneous kernels K(X,Y) are normally defined for
%     non-negative data only. VL_HOMKERMAP defines them for both
%     positive and negative data by using the definition
%     SIGN(X)SIGN(Y)K(ABS(X),ABS(Y)) -- note that other extensions are
%     possible as well (see [2]).
%
%   REFERENCES::
%     [1] A. Vedaldi and A. Zisserman
%     `Efficient Additive Kernels via Explicit Feature Maps',
%     Proc. CVPR, 2010.
%
%     [2] A. Vedaldi and A. Zisserman
%     `Efficient Additive Kernels via Explicit Feature Maps',
%     PAMI, 2011 (submitted).
%
%   See also: VL_HELP().

% Authors: Andrea Vedaldi

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: yihui-he

repo: panorama

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