text:
%% DEMO_febio_0055_clot_tube_slide_pressure.m
% Below is a demonstration for:
%
% * Building geometry for a spherical blob with tetrahedral elements
% which is being aspirated into a tube.
% This demo consists off:
% * Defining the boundary conditions
% * Coding the febio structure
% * Running the model
% * Importing and visualizing the displacement results

%% Keywords
%
% * febio_spec version 3.0
% * febio, FEBio
% * blood clot
% * contact, sliding, friction
% * rigid body constraints
% * hexahedral elements, hex8
% * quadrilaterl elements, quad4
% * shell elements
% * sphere
% * hyperelastic, Ogden
% * displacement logfile
% * stress logfile

%%

clear; close all; clc;

%% Plot settings
fontSize=15;
faceAlpha1=0.8;
faceAlpha2=0.3;
markerSize=40;
lineWidth=3;
cMap=[1 0.5 0.4; 0.9 0.3 0.27; 0.8 0.2 0.18; 0.7 0.1 0.09; 0.6 0 0; 0.5 0 0; 0.4 0 0;];
[cMap]=resampleColormap(cMap,250);

%% Control parameters

% Path names
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');

% Defining file names
febioFebFileNamePart='tempModel';
febioFebFileName=fullfile(savePath,[febioFebFileNamePart,'.feb']); %FEB file name
febioLogFileName=fullfile(savePath,[febioFebFileNamePart,'.txt']); %FEBio log file name
febioLogFileName_disp=[febioFebFileNamePart,'_disp_out.txt']; %Log file name for exporting displacement
febioLogFileName_strainEnergy=[febioFebFileNamePart,'_energy_out.txt']; %Log file name for exporting strain energy density

% Sphere parameters
sphereRadius=3;%
numElementsMantel=6;

% Ground plate parameters
tubeRadius=sphereRadius.*[1 0.1];
tubeAngle=3*(pi/180);
tubeLength=abs(diff(tubeRadius))/tan(tubeAngle);

% Material parameter set
c1=1e-4; %Shear-modulus-like parameter MPa
m1=2; %Material parameter setting degree of non-linearity
k_factor=10; %Bulk modulus factor
k=c1*k_factor; %Bulk modulus
d=1e-9; %Density

% FEA control settings
max_refs=25; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=15; %Optimum number of iterations
max_retries=25; %Maximum number of retires
symmetric_stiffness=0;
min_residual=1e-20;

timeTotal=1; %Analysis time
numTimeSteps=100; %Number of time steps desired
step_size=timeTotal/numTimeSteps;
dtmin=(timeTotal/numTimeSteps)/100; %Minimum time step size
dtmax=(timeTotal/numTimeSteps)*10; %Maximum time step size
analysisType='DYNAMIC';

%Contact parameters
contactPenalty=10;
laugon=0;
minaug=1;
maxaug=10;
fric_coeff=0.1;

%Specifying load
pressureValue=7.5e-5;

%% Creating model geometry and mesh
%

%Control settings
cPar.sphereRadius=sphereRadius;
cPar.coreRadius=cPar.sphereRadius/2;
cPar.numElementsMantel=numElementsMantel;
cPar.numElementsCore=round(numElementsMantel*1.5);
cPar.outputStructType=2;
cPar.makeHollow=0;
cPar.cParSmooth.n=25;

%Creating sphere
[meshOutput]=hexMeshSphere(cPar);

% Access model element and patch data
Fb_blob=meshOutput.facesBoundary;
Cb_blob=meshOutput.boundaryMarker;
V_blob=meshOutput.nodes;
E_blob=meshOutput.elements;

%%
% Visualize blob mesh

hFig=cFigure;
subplot(1,2,1); hold on;
gpatch(Fb_blob,V_blob,Cb_blob,'k',0.8);
patchNormPlot(Fb_blob,V_blob);
axisGeom(gca,fontSize);
colormap(gjet); icolorbar;
camlight headlight;

hs=subplot(1,2,2); hold on;
title('Cut view of solid mesh','FontSize',fontSize);
optionStruct.hFig=[hFig hs];
gpatch(Fb_blob,V_blob,'kw','none',0.25);
meshView(meshOutput,optionStruct);
axisGeom(gca,fontSize);
drawnow;

%% Creating tube model
%

pointSpacingBlob=mean(patchEdgeLengths(Fb_blob,V_blob));
pointSpacingTube=pointSpacingBlob;

rEnd=sphereRadius+(sphereRadius.*((sphereRadius-tubeRadius(2))/tubeLength));
V_curve_tube=[sphereRadius rEnd 0; -tubeLength tubeRadius(2) 0;];

nResample=ceil(max(pathLength(V_curve_tube))./pointSpacingTube);
V_curve_tube=evenlySampleCurve(V_curve_tube,nResample,'pchip',0);

cPar.closeLoopOpt=1;
cPar.numSteps=[]; %If empty the number of steps is derived from point spacing of input curve
cPar.w=[1 0 0];
[F_tube,V_tube]=polyRevolve(V_curve_tube,cPar);
[F_tube,V_tube]=mergeVertices(F_tube,V_tube);
center_of_mass_tube=mean(V_tube,1);

%% Join model node sets

V=[V_blob; V_tube; ];
F_tube=F_tube+size(V_blob,1);

%%
% Visualizing model

cFigure; hold on;
gtitle('Model components',fontSize);
hl(1)=gpatch(Fb_blob,V,'rw','k',0.8);
hl(2)=gpatch(F_tube,V,'kw','k',0.5);
legend(hl,{'Blob','Tube'}); clear hl;
axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Get contact surfaces
%

F_contact_blob=Fb_blob;

%%
% Visualize contact surfaces

cFigure; hold on;
title('Tube blob contact pair','fontsize',fontSize);
hl(1)=gpatch(F_tube,V,'rw','k',0.8);
patchNormPlot(F_tube,V);
hl(2)=gpatch(F_contact_blob,V,'kw','k',0.5);
patchNormPlot(F_contact_blob,V);
legend(hl,{'Master','Slave'}); clear hl;
axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Get pressure surface

N=patchNormal(F_contact_blob,V);
x=[1 0 0];
D=dot(N,x(ones(size(N,1),1),:),2);

logicFace=D>-1e-6; %Logic for current face set
F_pressure=F_contact_blob(logicFace,:); %The current face set

%%
% Visualizing boundary conditions. Markers plotted on the semi-transparent
% model denote the nodes in the various boundary condition lists.

cFigure; hold on;
title('Pressure surface','fontsize',fontSize);

gpatch(F_pressure,V,'kw','k',0.5);
patchNormPlot(F_pressure,V);

% legend(hl,{'Master','Slave'}); clear hl;
axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Defining the FEBio input structure
% See also |febioStructTemplate| and |febioStruct2xml| and the FEBio user
% manual.

%Get a template with default settings
[febio_spec]=febioStructTemplate;

%febio_spec version
febio_spec.ATTR.version='3.0';

%Module section
febio_spec.Module.ATTR.type='solid';

%Control section
febio_spec.Control.analysis=analysisType;
febio_spec.Control.time_steps=numTimeSteps;
febio_spec.Control.step_size=timeTotal/numTimeSteps;
febio_spec.Control.solver.max_refs=max_refs;
febio_spec.Control.solver.max_ups=max_ups;
febio_spec.Control.solver.symmetric_stiffness=symmetric_stiffness;
febio_spec.Control.time_stepper.dtmin=dtmin;
febio_spec.Control.time_stepper.dtmax=dtmax;
febio_spec.Control.time_stepper.max_retries=max_retries;
febio_spec.Control.time_stepper.opt_iter=opt_iter;

%Material section
materialName1='Material1';
febio_spec.Material.material{1}.ATTR.name=materialName1;
febio_spec.Material.material{1}.ATTR.type='Ogden unconstrained';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.c1=c1;
febio_spec.Material.material{1}.m1=m1;
febio_spec.Material.material{1}.c2=c1;
febio_spec.Material.material{1}.m2=-m1;
febio_spec.Material.material{1}.cp=k;
febio_spec.Material.material{1}.density=d;

materialName2='Material2';
febio_spec.Material.material{2}.ATTR.name=materialName2;
febio_spec.Material.material{2}.ATTR.type='rigid body';
febio_spec.Material.material{2}.ATTR.id=2;
febio_spec.Material.material{2}.density=1;
febio_spec.Material.material{2}.center_of_mass=center_of_mass_tube;

%Mesh section
% -> Nodes
febio_spec.Mesh.Nodes{1}.ATTR.name='nodeSet_all'; %The node set name
febio_spec.Mesh.Nodes{1}.node.ATTR.id=(1:size(V,1))'; %The node id's
febio_spec.Mesh.Nodes{1}.node.VAL=V; %The nodel coordinates

% -> Elements
partName1='Part1';
febio_spec.Mesh.Elements{1}.ATTR.name=partName1; %Name of this part
febio_spec.Mesh.Elements{1}.ATTR.type='hex8'; %Element type
febio_spec.Mesh.Elements{1}.elem.ATTR.id=(1:1:size(E_blob,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=E_blob; %The element matrix

partName2='Part2';
febio_spec.Mesh.Elements{2}.ATTR.name=partName2; %Name of this part
febio_spec.Mesh.Elements{2}.ATTR.type='quad4'; %Element type
febio_spec.Mesh.Elements{2}.elem.ATTR.id=size(E_blob,1)+(1:1:size(F_tube,1))'; %Element id's
febio_spec.Mesh.Elements{2}.elem.VAL=F_tube; %The element matrix

% -> Surfaces
surfaceName1='Surface1_contact';
febio_spec.Mesh.Surface{1}.ATTR.name=surfaceName1;
febio_spec.Mesh.Surface{1}.quad4.ATTR.id=(1:1:size(F_tube,1))';
febio_spec.Mesh.Surface{1}.quad4.VAL=F_tube;

surfaceName2='Surface2_contact';
febio_spec.Mesh.Surface{2}.ATTR.name=surfaceName2;
febio_spec.Mesh.Surface{2}.quad4.ATTR.id=(1:1:size(F_contact_blob,1))';
febio_spec.Mesh.Surface{2}.quad4.VAL=F_contact_blob;

surfaceName3='Surface3_pressure';
febio_spec.Mesh.Surface{3}.ATTR.name=surfaceName3;
febio_spec.Mesh.Surface{3}.quad4.ATTR.id=(1:1:size(F_pressure,1))';
febio_spec.Mesh.Surface{3}.quad4.VAL=F_pressure;

% -> Surface pairs
contactPairName='Contact1';
febio_spec.Mesh.SurfacePair{1}.ATTR.name=contactPairName;
febio_spec.Mesh.SurfacePair{1}.primary=surfaceName2;
febio_spec.Mesh.SurfacePair{1}.secondary=surfaceName1;

%MeshDomains section
febio_spec.MeshDomains.SolidDomain.ATTR.name=partName1;
febio_spec.MeshDomains.SolidDomain.ATTR.mat=materialName1;

febio_spec.MeshDomains.ShellDomain.ATTR.name=partName2;
febio_spec.MeshDomains.ShellDomain.ATTR.mat=materialName2;

%Loads
febio_spec.Loads.surface_load{1}.ATTR.type='pressure';
febio_spec.Loads.surface_load{1}.ATTR.surface=surfaceName3;
febio_spec.Loads.surface_load{1}.pressure.ATTR.lc=1;
febio_spec.Loads.surface_load{1}.pressure.VAL=pressureValue;
febio_spec.Loads.surface_load{1}.symmetric_stiffness=1;

%Rigid section
% -> Prescribed rigid body boundary conditions
febio_spec.Rigid.rigid_constraint{1}.ATTR.name='RigidFix_1';
febio_spec.Rigid.rigid_constraint{1}.ATTR.type='fix';
febio_spec.Rigid.rigid_constraint{1}.rb=2;
febio_spec.Rigid.rigid_constraint{1}.dofs='Rx,Ry,Rz,Ru,Rv,Rw';

%Contact section
febio_spec.Contact.contact{1}.ATTR.type='sliding-elastic';
febio_spec.Contact.contact{1}.ATTR.surface_pair=contactPairName;
febio_spec.Contact.contact{1}.two_pass=0;
febio_spec.Contact.contact{1}.laugon=laugon;
febio_spec.Contact.contact{1}.tolerance=0.2;
febio_spec.Contact.contact{1}.gaptol=0;
febio_spec.Contact.contact{1}.minaug=minaug;
febio_spec.Contact.contact{1}.maxaug=maxaug;
febio_spec.Contact.contact{1}.search_tol=0.01;
febio_spec.Contact.contact{1}.search_radius=0.01*sqrt(sum((max(V,[],1)-min(V,[],1)).^2,2));
febio_spec.Contact.contact{1}.symmetric_stiffness=0;
febio_spec.Contact.contact{1}.auto_penalty=1;
febio_spec.Contact.contact{1}.penalty=contactPenalty;
febio_spec.Contact.contact{1}.fric_coeff=fric_coeff;

%LoadData section
% -> load_controller
febio_spec.LoadData.load_controller{1}.ATTR.id=1;
febio_spec.LoadData.load_controller{1}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{1}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{1}.points.point.VAL=[0 0; 1 1];

%Output section
% -> log file
febio_spec.Output.logfile.ATTR.file=febioLogFileName;
febio_spec.Output.logfile.node_data{1}.ATTR.file=febioLogFileName_disp;
febio_spec.Output.logfile.node_data{1}.ATTR.data='ux;uy;uz';
febio_spec.Output.logfile.node_data{1}.ATTR.delim=',';

febio_spec.Output.logfile.element_data{1}.ATTR.file=febioLogFileName_strainEnergy;
febio_spec.Output.logfile.element_data{1}.ATTR.data='sed';
febio_spec.Output.logfile.element_data{1}.ATTR.delim=',';
febio_spec.Output.logfile.element_data{1}.VAL=1:size(E_blob,1);

%% Quick viewing of the FEBio input file structure
% The |febView| function can be used to view the xml structure in a MATLAB
% figure window.

%%
% |febView(febio_spec); %Viewing the febio file|

%% Exporting the FEBio input file
% Exporting the febio_spec structure to an FEBio input file is done using
% the |febioStruct2xml| function.

febioStruct2xml(febio_spec,febioFebFileName); %Exporting to file and domNode

%% Running the FEBio analysis
% To run the analysis defined by the created FEBio input file the
% |runMonitorFEBio| function is used. The input for this function is a
% structure defining job settings e.g. the FEBio input file name. The
% optional output runFlag informs the user if the analysis was run
% succesfully.

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.runMode='external';%'internal';

[runFlag]=runMonitorFEBio(febioAnalysis);%START FEBio NOW!!!!!!!!

%% Import FEBio results

if 1%runFlag==1 %i.e. a succesful run

    % Importing nodal displacements from a log file
    [time_mat, N_disp_mat,~]=importFEBio_logfile(fullfile(savePath,febioLogFileName_disp)); %Nodal displacements
    time_mat=[0; time_mat(:)]; %Time

    N_disp_mat=N_disp_mat(:,2:end,:);
    sizImport=size(N_disp_mat);
    sizImport(3)=sizImport(3)+1;
    N_disp_mat_n=zeros(sizImport);
    N_disp_mat_n(:,:,2:end)=N_disp_mat;
    N_disp_mat=N_disp_mat_n;
    DN=N_disp_mat(:,:,end);
    DN_magnitude=sqrt(sum(DN.^2,2));
    V_def=V+DN;
    V_DEF=N_disp_mat+repmat(V,[1 1 size(N_disp_mat,3)]);
    X_DEF=V_DEF(:,1,:);
    Y_DEF=V_DEF(:,2,:);
    Z_DEF=V_DEF(:,3,:);
    %     [CF]=vertexToFaceMeasure(Fb_all,DN_magnitude);

    %%
    % Importing element strain energies from a log file
    [~,E_energy,~]=importFEBio_logfile(fullfile(savePath,febioLogFileName_strainEnergy)); %Element stresses

    %Remove nodal index column
    E_energy=E_energy(:,2:end,:);

    %Add initial state i.e. zero displacement
    sizImport=size(E_energy);
    sizImport(3)=sizImport(3)+1;
    E_energy_mat_n=zeros(sizImport);
    E_energy_mat_n(:,:,2:end)=E_energy;
    E_energy=E_energy_mat_n;

    %%

    indBlob=unique(Fb_blob(:));
    t=linspace(0,2*pi,250)';

    V_def_blob=V(indBlob,:)+N_disp_mat(indBlob,:,end);

    [~,indMin]=min(V_def_blob(:,1));
    [~,indMax]=max(V_def_blob(:,1));

    xEnd=V_def_blob(indMin,1);
    xStart=V_def_blob(indMax,1);
    rEnd=sphereRadius+(xEnd.*((sphereRadius-tubeRadius(2))/tubeLength));
    rStart=sphereRadius+(xStart.*((sphereRadius-tubeRadius(2))/tubeLength));

    xMid=mean([xStart xEnd]);%sum([rStart rEnd].*[xStart xEnd])./sum([rStart rEnd]);
    rMid=sphereRadius+(xMid.*((sphereRadius-tubeRadius(2))/tubeLength));

    V_plot_xEnd=[xEnd*ones(size(t)) rEnd*cos(t) rEnd*sin(t)];
    V_plot_xMid=[xMid*ones(size(t)) rMid*cos(t) rMid*sin(t)];
    V_plot_xStart=[xStart*ones(size(t)) rStart*cos(t) rStart*sin(t)];

   %%
    % Plotting the simulated results using |anim8| to visualize and animate
    % deformations

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; hold on;
    ht=gtitle(['Radial stretch: ',num2str(rMid/sphereRadius)]);
    hp1=gpatch(Fb_blob,V_def,DN_magnitude,'none',1); %Add graphics object to animate

    hp2=plotV(V_plot_xEnd  ,'r-','LineWidth',3);
    hp3=plotV(V_plot_xMid  ,'r-','LineWidth',3);
    hp4=plotV(V_plot_xStart,'r-','LineWidth',3);

    gpatch(F_tube,V_def,'kw','none',0.25); %Add graphics object to animate
    axisGeom(gca,fontSize);
    colormap(cMap); colorbar;
    caxis([0 max(DN_magnitude(:))]); caxis manual;
    axis([min(X_DEF(:)) max(X_DEF(:)) min(Y_DEF(:)) max(Y_DEF(:)) min(Z_DEF(:)) max(Z_DEF(:))]);
    camlight headlight; lighting gouraud;
    view(0,0);
%     view(-30,30); zoom(1.5);
    axis off;
    drawnow;

    LMid=1;

    % Set up animation features
    animStruct.Time=time_mat; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments
        DN=N_disp_mat(:,:,qt); %Current displacement
        DN_magnitude=sqrt(sum(DN.^2,2));
        V_def=V+DN; %Current nodal coordinates

        V_def_blob=V_def(indBlob,:);
        [~,indMin]=min(V_def_blob(:,1));
        [~,indMax]=max(V_def_blob(:,1));

        xEnd=V_def_blob(indMin,1);
        xStart=V_def_blob(indMax,1);
        rEnd=sphereRadius+(xEnd.*((sphereRadius-tubeRadius(2))/tubeLength));
        rStart=sphereRadius+(xStart.*((sphereRadius-tubeRadius(2))/tubeLength));

        xMid=mean([xStart xEnd]);%sum([rStart rEnd].*[xStart xEnd])./sum([rStart rEnd]);
        rMid=sphereRadius+(xMid.*((sphereRadius-tubeRadius(2))/tubeLength));

        LMid=min(LMid,rMid/sphereRadius);

        V_plot_xEnd=[xEnd*ones(size(t)) rEnd*cos(t) rEnd*sin(t)];
        V_plot_xMid=[xMid*ones(size(t)) rMid*cos(t) rMid*sin(t)];
        V_plot_xStart=[xStart*ones(size(t)) rStart*cos(t) rStart*sin(t)];

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp1 hp1 hp2 hp2 hp2 hp3 hp3 hp3 hp4 hp4 hp4 ht]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData','XData','YData','ZData','XData','YData','ZData','XData','YData','ZData','String'}; %Properties of objects to animate
        animStruct.Set{qt}={V_def,DN_magnitude,...
            V_plot_xEnd(:,1),V_plot_xEnd(:,2),V_plot_xEnd(:,3),...
            V_plot_xMid(:,1),V_plot_xMid(:,2),V_plot_xMid(:,3),...
            V_plot_xStart(:,1),V_plot_xStart(:,2),V_plot_xStart(:,3),...
            ['Radial stretch: ',num2str(rMid/sphereRadius)]}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature

    drawnow;

end

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*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


text:
%Walsh-Harmard Transform with Tow level Discrete Wavelete Transform for
% Color image Decompression
% Designed by  Mohammed M. Siddeq
% Data 2012-2-22
% Email :- mamadmmx76@yahoo.com
%
% this program is used for Decompress grayscale images by using :
% INPUT\  Header : this parameter contains all infmration about Color Compressed file
%
%OUTPUT\  Im : Decoded Color image from "Header"

clear;
Path_Name='C:\WWT\images\2.wwt';
X=load(Path_Name, '-mat'); % Read compress data from the file
Header=X.Header;
H1=Header(1).H1;
H2=Header(2).H2;
H3=Header(3).H3;
Y=Walsh_DWT_Decoding(H1);% Apply Decompression on each layer
Cb=Walsh_DWT_Decoding(H2);
Cr=Walsh_DWT_Decoding(H3);

%% Collect all layers in one matrix Ycbcr
ycbcr(:,:,1)=Y(:,:);
ycbcr(:,:,2)=Cb(:,:);
ycbcr(:,:,3)=Cr(:,:);

Im= ycbcr2rgb(ycbcr); % Convert YCbCr format to RGB

imshow(uint8(Im));

author: Sable

repo: mcbench-benchmarks


text:
function obj1 = plus(obj1,obj2)

for i = 1:length(obj1.u)
  obj1.u{i} = obj1.u{i} + obj2.u{i};
end

end

author: mtex-toolbox

repo: mtex

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


text:
% Copyright (C) Daphne Koller, Stanford University, 2012

function [num mv] = randi(arg1,arg2,arg3)

num = -1;
persistent x_i;
persistent p1;
persistent p2;
if(isempty(x_i))
  x_i = 1;
  p1 = 160481183;
  p2 = 179424673;
end
mv=p2;
if(ischar(arg1)==1)
  if(strcmp(arg1,'seed'))
    if(nargin>1)
      x_i = arg2;
      num = 0;
    else
      x_i=1;
      num=0;
    end
  else
    'Unrecognized option. The only accepted option to this random library is -seed-.'
  end
else
  if(arg1>p2)
    'Max too high, range cutoff at 1 million'
  end
  if(nargin>1)
    if(nargin==2)
      arg3=arg2;
    end
    num = zeros(arg2,arg3);
    for i=1:arg2
      for j = 1:arg3
        x_i = mod(x_i*(p1+1)+p1,p2);
        num(i,j) = mod(x_i,arg1)+1;
      end
    end
  else
    x_i = mod(x_i*(p1+1)+p1,p2);
    num=mod(x_i,arg1)+1;
  end
end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function f=comp_idgt_fac(coef,gf,L,a,M)
%COMP_IDGT_FAC  Full-window factorization of a Gabor matrix.
%   Usage:  f=comp_idgt_fac(c,g,a,M)
%
%   Input parameters:
%         c     : M x N array of coefficients.
%         gf    : Factorization of window (from facgabm).
%         a     : Length of time shift.
%         M     : Number of frequency shifts.
%   Output parameters:
%         f     : Reconstructed signal.
%
%   Do not call this function directly, use IDGT.
%   This function does not check input parameters!
%
%   If input is a matrix, the transformation is applied to
%   each column.
%
%   This function does not handle multidimensional data, take care before
%   you call it.
%
%   References: so07-2 st98-8

%   AUTHOR : Peter L. Søndergaard.
%   TESTING: OK
%   REFERENCE: OK

% Calculate the parameters that was not specified.
N=L/a;
b=L/M;

R=prod(size(gf))/L;

W=prod(size(coef))/(M*N*R);

N=L/a;
b=L/M;

[c,h_a,h_m]=gcd(a,M);
h_a=-h_a;
p=a/c;
q=M/c;
d=N/q;

ff=zeros(p,q*W,c,d,assert_classname(coef,gf));
C=zeros(q*R,q*W,c,d,assert_classname(coef,gf));
f=zeros(L,W,assert_classname(coef,gf));

% Apply ifft to the coefficients.
%coef=ifft(reshape(coef,M,N*W))*sqrt(M);
coef=ifft(coef)*sqrt(M);

% Set up the small matrices

coef=reshape(coef,M,N,R,W);

if p==1

  for rw=0:R-1
    for w=0:W-1
      for s=0:d-1
	for l=0:q-1
	  for u=0:q-1
	    C(u+1+rw*q,l+1+w*q,:,s+1)=coef((1:c)+l*c,mod(u+s*q+l,N)+1,rw+1,w+1);
	  end;
	end;
      end;
    end;
  end;
else
  % Rational oversampling
  for rw=0:R-1
    for w=0:W-1
      for s=0:d-1
	for l=0:q-1
	  for u=0:q-1
	    C(u+1+rw*q,l+1+w*q,:,s+1)=coef((1:c)+l*c,mod(u+s*q-l*h_a,N)+1,rw+1,w+1);
	  end;
	end;
      end;
    end;
  end;
end;

% FFT them
if d>1
  C=fft(C,[],4);
end;

% Multiply them
for r=0:c-1
  for s=0:d-1
    CM=reshape(C(:,:,r+1,s+1),q*R,q*W);
    GM=reshape(gf(:,r+s*c+1),p,q*R);

    ff(:,:,r+1,s+1)=GM*CM;
  end;
end;

% Inverse FFT
if d>1
  ff=ifft(ff,[],4);
end;

% Place the result
if p==1

  for s=0:d-1
    for w=0:W-1
      for l=0:q-1
	f((1:c)+mod(s*M+l*a,L),w+1)=reshape(ff(1,l+1+w*q,:,s+1),c,1);
      end;
    end;
  end;

else
  % Rational oversampling
  for w=0:W-1
    for s=0:d-1
      for l=0:q-1
	for k=0:p-1
	  f((1:c)+mod(k*M+s*p*M-l*h_a*a,L),w+1)=reshape(ff(k+1,l+1+w*q,:,s+1),c,1);
	end;
      end;
    end;
  end;

end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function o = npRegGet(options,name,default,flag)
%NPREGGET Get NPREG OPTIONS parameters.
%   VAL = NPREGGET(OPTIONS,'NAME') extracts the value of the named
%   parameter from npReg options structure OPTIONS, returning an empty
%   matrix if the parameter value is not specified in OPTIONS.  It is
%   sufficient to type only the leading characters that uniquely identify
%   the parameter.  Case is ignored for parameter names.  [] is a valid
%   OPTIONS argument.
%
%   VAL = NPREGGET(OPTIONS,'NAME',DEFAULT) extracts the named parameter
%   as above, but returns DEFAULT if the named parameter is not specified
%   (is []) in OPTIONS.  For example
%
%     val = npRegGet(opts,'UDiffTol',1e-4);
%
%   returns val = 1e-4 if the UDiffTol property is not specified in opts.
%
%   See also NPREGSET.
%
% author: Nathan D. Cahill
% email: nathan.cahill@rit.edu
% affiliation: Rochester Institute of Technology
% date: January 2014
% licence: GNU GPL v3 licence.
%
% Copyright Nathan D. Cahill
% Code available from https://github.com/tomdoel/npReg

% undocumented usage for fast access with no error checking
if (nargin == 4) && isequal(flag,'fast')
    o = npRegGetFast(options,name,default);
    return
end

if nargin < 2
    error('MATLAB:npRegGet:NotEnoughInputs', 'Not enough input arguments.');
end
if nargin < 3
    default = [];
end

if ~isempty(options) && ~isa(options,'struct')
    error('MATLAB:npRegGet:Arg1NotStruct',...
        'First argument must be an options structure created with NPREGSET.');
end

if isempty(options)
    o = default;
    return;
end

allfields = {'Display'};
try
    npRegFields = npRegOptionGetFields;
    allfields = [allfields; npRegFields];
catch
    lasterrstruct = lasterror;
    if strcmp(lasterrstruct.identifier, 'MATLAB:UndefinedFunction')
        % Function NPREGOPTIONGETFIELDS not found, so we assume npReg Toolbox not on path
        %   and use the "MATLAB only" struct.
        % clean up last error
        lasterr('');
    else
        rethrow(lasterror);
    end
end

Names = allfields;

name = deblank(name(:)'); % force this to be a row vector
j = find(strncmpi(name,Names,length(name)));
if isempty(j)               % if no matches
    error('MATLAB:npRegGet:InvalidPropName',...
        ['Unrecognized property name ''%s''.  ' ...
        'See NPREGSET for possibilities.'], name);
elseif length(j) > 1            % if more than one match
    % Check for any exact matches (in case any names are subsets of others)
    k = find(strcmpi(name,Names));
    if length(k) == 1
        j = k;
    else
        msg = sprintf('Ambiguous property name ''%s'' ', name);
        msg = [msg '(' Names{j(1),:}];
        for k = j(2:length(j))'
            msg = [msg ', ' Names{k,:}];
        end
        msg = sprintf('%s).', msg);
        error('MATLAB:npRegGet:AmbiguousPropName', msg);
    end
end

if any(strcmp(Names,Names{j,:}))
    o = options.(Names{j,:});
    if isempty(o)
        o = default;
    end
else
    o = default;
end

%------------------------------------------------------------------
function value = npRegGetFast(options,name,defaultopt)
%NPREGGETFAST Get NPREG OPTIONS parameter with no error checking.
%   VAL = NPREGGETFAST(OPTIONS,FIELDNAME,DEFAULTOPTIONS) will get the
%   value of the FIELDNAME from OPTIONS with no error checking or
%   fieldname completion. If the value is [], it gets the value of the
%   FIELDNAME from DEFAULTOPTIONS, another OPTIONS structure which is
%   probably a subset of the options in OPTIONS.
%

if isempty(options)
     value = defaultopt.(name);
     return;
end
% We need to know if name is a valid field of options, but it is faster to use
% a try-catch than to test if the field exists and if the field name is
% correct. If the options structure is from an older version of the
% toolbox, it could be missing a newer field.
try
    value = options.(name);
catch
    value = [];
    lasterr('');  % clean up last error
end

if isempty(value)
    value = defaultopt.(name);
end

author: tomdoel

repo: pulmonarytoolkit


text:
function installBayesFactor
% Install the bayesFactor class in Matlab.

lmm = which('linearmixedmodel');
if isempty(lmm)
    disp('Hmm. You don''t seem to have the the Statistics and Machine Learning Toolbox. This is not likely to work');
end

filename = mfilename('fullpath');
pth = fileparts(filename);
addpath(pth)
savepath;

% Store the default options as prefs on this machine.
% The user can change these with bf.options
bf.options;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function means = gtmlmean(net, data)
%GTMLMEAN Mean responsibility for data in a GTM.
%
%	Description
%	 MEANS = GTMLMEAN(NET, DATA) takes a GTM structure NET, and computes
%	the means of the responsibility  distributions for each data point in
%	DATA.
%
%	See also
%	GTM, GTMPOST, GTMLMODE
%

%	Copyright (c) Ian T Nabney (1996-2001)

% Check for consistency
errstring = consist(net, 'gtm', data);
if ~isempty(errstring)
  error(errstring);
end

R = gtmpost(net, data);
means = R*net.X;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function overlap = calcRectInt(A,B)
%
%each row is a rectangle.
% A(i,:) = [x y w h]
% B(j,:) = [x y w h]
% overlap(i,j) = area of intersection
% normoverlap(i,j) = overlap(i,j) / (area(i)+area(j)-overlap)
%
% Same as built-in rectint, but faster and uses less memory (since avoids repmat).

lena = size(A, 1);
lenb = size(B, 1);
lena = min(lena, lenb);
A = A(1:lena, :);
B = B(1:lena, :);
leftA = A(:,1);
bottomA = A(:,2);
rightA = leftA + A(:,3) - 1;
topA = bottomA + A(:,4) - 1;

leftB = B(:,1);
bottomB = B(:,2);
rightB = leftB + B(:,3) - 1;
topB = bottomB + B(:,4) - 1;

tmp = (max(0, min(rightA, rightB) - max(leftA, leftB)+1 )) .* (max(0, min(topA, topB) - max(bottomA, bottomB)+1 ));
    areaA = A(:,3) .* A(:,4);
    areaB = B(:,3) .* B(:,4);
    overlap = tmp./(areaA+areaB-tmp);
% if tmp > 0
%
%     overlap = tmp;
%
%     areaA = A(3) .* A(4);
%     areaB = B(3) .* B(4);
%     overlap = tmp./(areaA+areaB-tmp);
% else
%     overlap = 0;
% end

repo: STRCF


text:
function plin = linearize(p)
% LINEARIZE Linearize SDPVAR object
%
% h = LINEARIZE(p)
%
% Returns linearization p(double(x)) + dp(double(x))*(x-double(x))
% where x is the SDPVAR variables defining the polynomial p(x)
%
% See also SDPVAR, JACOBIAN

if isa(p,'double')
    plin = zeros(size(p));
    return
end

if is(p,'linear') & ~is(p,'compound')
    plin = p;
    return
end

x = recover(depends(p));
x0 = double(x);
p0 = double(p);

n = size(p,1);
m = size(p,2);

if min(n,m)>1
    plin = [];
    for i = 1:m
        plin = [plin p0(:,i)+double(jacobian(p(:,i),x))*(x-x0)];
    end
else
    plin = p0+double(jacobian(p,x))*(x-x0);
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function wcoef = FWT_PO(x,L,qmf)
% FWT_PO -- Forward Wavelet Transform (periodized, orthogonal)
%  Usage
%    wc = FWT_PO(x,L,qmf)
%  Inputs
%    x    1-d signal; length(x) = 2^J
%    L    Coarsest Level of V_0;  L << J
%    qmf  quadrature mirror filter (orthonormal)
%  Outputs
%    wc    1-d wavelet transform of x.
%
%  Description
%    1. qmf filter may be obtained from MakeONFilter
%    2. usually, length(qmf) < 2^(L+1)
%    3. To reconstruct use IWT_PO
%
%  See Also
%    IWT_PO, MakeONFilter
%
  [n,J] = dyadlength(x) ;
  wcoef = zeros(1,n) ;
  beta = ShapeAsRow(x);  %take samples at finest scale as beta-coeffts
  for j=J-1:-1:L
       alfa = DownDyadHi(beta,qmf);
       wcoef(dyad(j)) = alfa;
       beta = DownDyadLo(beta,qmf) ;
  end
  wcoef(1:(2^L)) = beta;
  wcoef = ShapeLike(wcoef,x);

%
% Copyright (c) 1993. Iain M. Johnstone
%

%
%  Part of Wavelab Version 850
%  Built Tue Jan  3 13:20:40 EST 2006
%  This is Copyrighted Material
%  For Copying permissions see COPYING.m
%  Comments? e-mail wavelab@stat.stanford.edu

author: BUNPC


text:
function linplus_test445 ( )

%*****************************************************************************80
%
%% TEST445 tests R8PBL_INDICATOR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  n = 9;
  mu = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST445\n' );
  fprintf ( 1, '  R8PBL_INDICATOR sets up a R8PBL indicator matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
  fprintf ( 1, '  Bandwidth MU = %d\n', mu );

  a = r8pbl_indicator ( n, mu );

  r8pbl_print ( n, mu, a, '  The R8PBL indicator matrix:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function planC = annotatePlanCForXNAT(planC, xhost,xexp,xproj,xsubj)
% function planC = annotatePlanCForXNAT(planC, xhost,xexp,xproj,xsubj)
%
% This function adds XNAT addressing metadata to an existing planC

if nargin < 3 || isempty(xproj)
    xproj = '';
    if isempty(xsubj)
        xsubj = '';
    end
end

indexS = planC{end};
planC{indexS.header}.xnatInfo = struct('host',xhost,'experiment_id',xexp,'project_id',xproj,'subject_id',xsubj);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ r, code, w, m ] = s_to_format ( s )

%*****************************************************************************80
%
%% S_TO_FORMAT reads a FORTRAN format from a string.
%
%  Discussion:
%
%    This routine will read as many characters as possible until it reaches
%    the end of the string, or encounters a character which cannot be
%    part of the format.  This routine is limited in its ability to
%    recognize FORTRAN formats.  In particular, we are only expecting
%    a single format specification, and cannot handle extra features
%    such as 'ES' and 'EN' codes, '5X' spacing, and so on.
%
%    Legal input is:
%
%       0 nothing
%       1 blanks
%       2 optional '('
%       3 blanks
%       4 optional repeat factor R
%       5 blanks
%       6 CODE ( 'A', 'B', 'E', 'F', 'G', 'I', 'L', 'O', 'Z', '*' )
%       7 blanks
%       8 width W
%       9 optional decimal point
%      10 optional mantissa M
%      11 blanks
%      12 optional ')'
%      13 blanks
%
%  Example:
%
%    S                 R   CODE   W    M
%
%    'I12              1   I      12   0
%    'E8.0'            1   E       8   0
%    'F10.5'           1   F      10   5
%    '2G14.6'          2   G      14   6
%    '*'               1   *      -1  -1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string containing the
%    data to be read.  Reading will begin at position 1 and
%    terminate at the end of the string, or when no more
%    characters can be read.
%
%    Output, integer R, the repetition factor, which defaults to 1.
%
%    Output, character CODE, the format code.
%
%    Output, integer W, the field width.
%
%    Output, integer M, the mantissa width.
%
  LEFT = 1;
  RIGHT = -1;

  state = 0;
  paren_sum = 0;
  pos = 0;
  s_length = s_len_trim ( s );

  r = 0;
  w = 0;
  code = '?';
  m = 0;

  while ( pos < s_length )

    pos = pos + 1;
    c = s(pos);
%
%  BLANK character:
%
    if ( c == ' ' )

      if ( state == 4 )
        state = 5;
      elseif ( state == 6 )
        state = 7;
      elseif ( state == 10 )
        state = 11;
      elseif ( state == 12 )
        state = 13;
      end
%
%  LEFT PAREN
%
    elseif ( c == '(' )

      if ( state < 2 )
        paren_sum = paren_sum + LEFT;
      else
        state = -1;
        break;
      end
%
%  DIGIT (R, F, or W)
%
    elseif ( ch_is_digit ( c ) )

      if ( state <= 3 )
        state = 4;
        r = ch_to_digit ( c );
      elseif ( state == 4 )
        d = ch_to_digit ( c );
        r = 10 * r + d;
      elseif ( state == 6 || state == 7 )
        if ( code == '*' )
          state = -1;
          break;
        end
        state = 8;
        w = ch_to_digit ( c );
      elseif ( state == 8 )
        d = ch_to_digit ( c );
        w = 10 * w + d;
      elseif ( state == 9 )
        state = 10;
        m = ch_to_digit ( c );
      elseif ( state == 10 )
        d = ch_to_digit ( c );
        m = 10 * m + d;
      else
        state = -1;
        break;
      end
%
%  DECIMAL POINT
%
    elseif ( c == '.' )

      if ( state == 8 )
        state = 9;
      else
        state = -1;
        break;
      end
%
%  RIGHT PAREN
%
    elseif ( c == ')' )

      paren_sum = paren_sum + RIGHT;

      if ( paren_sum ~= 0 )
        state = -1;
        break;
      end

      if ( state == 6 & code == '*' )
        state = 12;
      elseif ( 6 <= state )
        state = 12;
      else
        state = -1;
        break;
      end
%
%  Code
%
    elseif ( ch_is_format_code ( c ) )

      if ( state < 6 )
        state = 6;
        code = c;
      else
        state = -1;
        break;
      end
%
%  Unexpected character
%
    else

      state = -1;
      break;

    end

  end

  if ( paren_sum ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'S_TO_FORMAT - Fatal error!\n' );
    fprintf ( 1, '  Parentheses mismatch.\n' );
    error ( 'S_TO_FORMAT - Fatal error!' );
  end

  if ( state < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'S_TO_FORMAT - Fatal error!\n' );
    fprintf ( 1, '  Parsing error.\n' );
    error ( 'S_TO_FORMAT - Fatal error!' );
  end

  if ( r == 0 )
    r = 1;
  end

  return
end


text:
function EditrunGui
% Modifiable runGUI file
clc;
clear all;
fileName = 'zoom_strips_plot.mat';    %USER - ENTER FILENAME
fileData=load(fileName);
temp=fileData(1).temp;

f = figure('Visible','on',...
'Units','normalized',...
'Position',[0,0,1,1],...
'MenuBar','none',...
'NumberTitle','off');

%SENSE COMPUTER AND SET FILE DELIMITER
switch(computer)
    case 'MACI64',		char= '/';
    case 'GLNX86',  char='/';
    case 'PCWIN',	char= '\';
    case 'PCWIN64', char='\';
    case 'GLNXA64', char='/';
end

% start_path='c:\data\matlab_gui\speech_files';

% find speech files directory by going up one level and down one level
% on the directory chain; as follows:
    dir_cur=pwd; % this is the current Matlab exercise directory path
    s=regexp(dir_cur,'\'); % find the last '\' for the current directory
    s1=s(length(s)); % find last '\' character; this marks upper level directory
    dir_fin=strcat(dir_cur(1:s1),'speech_files'); % create new directory
    start_path=dir_fin; % save new directory for speech files location

Callbacks_zoom_strips_plot_GUI25(f,temp,start_path);    %USER - ENTER PROPER CALLBACK FILE
%panelAndButtonEdit(f, temp);       % Easy access to Edit Mode

% Note comment PanelandBUttonCallbacks(f,temp) if panelAndButtonEdit is to
% be uncommented and used
end

% 2 Panels -- one for speech processing parameters, one for graphics
% 1 Graphics Panel -- strips waveform plot
% 1 Text Box
% 8 Buttons
%   #1 - pushbutton - Get Speech Directory (speech_files)
%   #2 - popupmenu - array of speech files from directory speech_files
%   #3 - editable button - starting sample of strips waveform plot
%   #4 - editable button - number of samples per line
%   #5 - popupmenu - samples or seconds for horizontal axis
%   #6 - pushbutton - Play Speech File
%   #7 - pushbutton - Plot Speech File
%   #8 - pushbutton - Close GUI


text:
function favar_excelrecord6fcn(identified, IRFperiods, favar, labels, favar_fevd_estimates, pref)

% create the cell that will be saved on excel
fevdcell={};

% build preliminary elements: space between the tables
horzspace=repmat({''},2,5*(identified+1));
vertspace=repmat({''},IRFperiods+3,1);

%subsample
informationvariablestrings=favar.informationvariablestrings(1,favar.plotX_index);

% loop over variables (vertical dimension)
for ii=1:favar.npltX
    tempcell={};
    % loop over shocks (horizontal dimension)
    for jj=1:identified+1
        % create cell of fevd record for variable ii in response to shock jj
        temp=['part of ' informationvariablestrings{1,ii} ' fluctuation due to ' labels{jj,1} ' shocks'];
        fevd_ij=[temp {''} {''} {''};{''} {''} {''} {''};{''} {'lw. bound'} {'median'} {'up. bound'};[num2cell((1:IRFperiods)') num2cell((favar_fevd_estimates{ii,jj})')]];
        tempcell=[tempcell fevd_ij vertspace];
    end
    fevdcell=[fevdcell;horzspace;tempcell];
end
% trim
fevdcell=fevdcell(3:end,1:end-1);
% write in excel
warning off MATLAB:xlswrite:AddSheet
xlswrite(fullfile(pref.results_path, [pref.results_sub '.xlsx']),fevdcell,'favar_FEVD','B2');
warning on MATLAB:xlswrite:AddSheet

author: european-central-bank

repo: BEAR-toolbox


text:
function d = cols(X)

warning('This function is deprecated')

d = size(X,2);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [R,C,B] = circumradius(V,T)
  % CIRCUMRADIUS Return the circumradius of each triangle/tet element
  %
  % R = circumradius(V,T)
  % [R,C,B] = circumradius(V,T)
  %
  % Input:
  %   V  #V by dim list of vertex positions
  %   T  #T by simplex-size list of tet indices
  % Output:
  %   R  #T by 1 list of simplex circumradii
  %   C  #T by dim list of simplex circumcenters
  %   B  #T by simplex-size list of barycentric coordinates so that:
  %     C(t,:) = B(t,:) * V(T(t,:),:)
  %
  % Known issues:
  %   B output only supported for triangles
  %

  switch size(T,2)
  case 4
    d_14 = sqrt(sum((V(T(:,1),:)-V(T(:,4),:)).^2,2));
    d_23 = sqrt(sum((V(T(:,2),:)-V(T(:,3),:)).^2,2));
    d_24 = sqrt(sum((V(T(:,2),:)-V(T(:,4),:)).^2,2));
    d_34 = sqrt(sum((V(T(:,3),:)-V(T(:,4),:)).^2,2));
    d_12 = sqrt(sum((V(T(:,1),:)-V(T(:,2),:)).^2,2));
    d_13 = sqrt(sum((V(T(:,1),:)-V(T(:,3),:)).^2,2));
    R = sqrt((d_14.*d_23+d_13.*d_24-d_12.*d_34).*(d_14.*d_23-d_13.*d_24+d_12.*d_34).*(-d_14.*d_23+d_13.*d_24+d_12.*d_34).*(d_14.*d_23+d_13.*d_24+d_12.*d_34))./(sqrt(2)*sqrt(-2.*d_34.^2.*d_12.^4-2.*d_34.^4.*d_12.^2-2.*d_13.^2.*d_23.^2.*d_12.^2+2.*d_14.^2.*d_23.^2.*d_12.^2+2.*d_13.^2.*d_24.^2.*d_12.^2-2.*d_14.^2.*d_24.^2.*d_12.^2+2.*d_13.^2.*d_34.^2.*d_12.^2+2.*d_14.^2.*d_34.^2.*d_12.^2+2.*d_23.^2.*d_34.^2.*d_12.^2+2.*d_24.^2.*d_34.^2.*d_12.^2-2.*d_14.^2.*d_23.^4-2.*d_13.^2.*d_24.^4-2.*d_14.^4.*d_23.^2+2.*d_13.^2.*d_14.^2.*d_23.^2-2.*d_13.^4.*d_24.^2+2.*d_13.^2.*d_14.^2.*d_24.^2+2.*d_13.^2.*d_23.^2.*d_24.^2+2.*d_14.^2.*d_23.^2.*d_24.^2-2.*d_13.^2.*d_14.^2.*d_34.^2+2.*d_14.^2.*d_23.^2.*d_34.^2+2.*d_13.^2.*d_24.^2.*d_34.^2-2.*d_23.^2.*d_24.^2.*d_34.^2));
    assert(nargout == 1);
  case 3
    % http://www.mathworks.com/matlabcentral/fileexchange/17300-circumcircle-of-a-triangle/content/circumcircle.m

    %compute the length of sides (AB, BC and CA) of the triangle
    l = [ ...
      sqrt(sum((V(T(:,2),:)-V(T(:,3),:)).^2,2)) ...
      sqrt(sum((V(T(:,3),:)-V(T(:,1),:)).^2,2)) ...
      sqrt(sum((V(T(:,1),:)-V(T(:,2),:)).^2,2)) ...
      ];
    dblA = doublearea(V,T);
    %use formula: R=abc/(4*area) to compute the circum radius
    R = prod(l,2)./(2*dblA);
    %compute the barycentric coordinates of the circum center
    B= [ ...
      l(:,1).^2.*(-l(:,1).^2+l(:,2).^2+l(:,3).^2) ...
      l(:,2).^2.*( l(:,1).^2-l(:,2).^2+l(:,3).^2) ...
      l(:,3).^2.*( l(:,1).^2+l(:,2).^2-l(:,3).^2)];
    % normalize
    B = bsxfun(@rdivide,B,sum(B,2));
    %convert to the real coordinates
    C = zeros(size(T,1),size(V,2));
    for d = 1:size(V,2);
      C(:,d) = sum(B.*[V(T(:,1),d) V(T(:,2),d) V(T(:,3),d)],2);
    end

end

author: alecjacobson

repo: gptoolbox


text:
function [out1] = points_on_line(inp1, inp2, inp3, inp4, inp5)
out1 = [];
ns = getNamespace();
if isempty(ns)
    return;
end
if strcmp(ns, 'AtlasViewerGUI')
    if nargin == 0
        [out1] = points_on_line_AtlasViewerGUI();
    elseif nargin == 1
        [out1] = points_on_line_AtlasViewerGUI(inp1);
    elseif nargin == 2
        [out1] = points_on_line_AtlasViewerGUI(inp1, inp2);
    elseif nargin == 3
        [out1] = points_on_line_AtlasViewerGUI(inp1, inp2, inp3);
    elseif nargin == 4
        [out1] = points_on_line_AtlasViewerGUI(inp1, inp2, inp3, inp4);
    elseif nargin == 5
        [out1] = points_on_line_AtlasViewerGUI(inp1, inp2, inp3, inp4, inp5);
    end
elseif strcmp(ns, 'Homer3') || strcmp(ns, 'DataTreeClass')
    if nargin == 0
        [out1] = points_on_line_Homer3();
    elseif nargin == 1
        [out1] = points_on_line_Homer3(inp1);
    elseif nargin == 2
        [out1] = points_on_line_Homer3(inp1, inp2);
    elseif nargin == 3
        [out1] = points_on_line_Homer3(inp1, inp2, inp3);
    elseif nargin == 4
        [out1] = points_on_line_Homer3(inp1, inp2, inp3, inp4);
    end
end

% ---------------------------------------------------------
function [p3] = points_on_line_AtlasViewerGUI(p1, p2, crop, mode, gridsize)

%
% USAGE:
%
%    p3 = points_on_line_AtlasViewerGUI(p1,p2,crop,mode)
%
% DESCRIPTION:
%
%    Function takes 2 points in 3d, p1 and p2, and a crop length
%    and finds a point on the line formed by extending p1,p2 by the
%    amount (or fraction of p1-p2 distance) crop extends away from p1.
%
%    Positive crop crops p1,p2 by the crop amount.
%    Negative crop extends p1,p2 by the crop amount.
%
%    mode = 'relative' means crop is the distance p3 extends away from p1 relative
%    to the distance between p1 and p2. This is the default if no mode is provided.
%
%    mode = 'absolute' means crop is the absolute distance (in unit length) that
%    p3 extends away from p1.
%
%    mode = 'all' returns the set of all points at intervals 1/crop on the line
%    segment between p1 and p2.
%
% AUTHOR: Jay Dubb (jdubb@nmr.mgh.harvard.edu)
% DATE:   7/7/2012
%

if ~exist('mode','var') || isempty(mode)
    mode = 'relative';
end
if ~exist('gridsize','var') || isempty(gridsize)
    gridsize = 0;
end

p3 = p1;

if isempty(p1) || isempty(p2)
    return;
end
if all(p1==p2)
    return;
end

if ~exist('crop','var') || isempty(crop)
    crop = set_line_crop_AtlasViewerGUI(p2, p1, gridsize);
end

dx0 = p1(1)-p2(1);
dy0 = p1(2)-p2(2);
dz0 = p1(3)-p2(3);
dr = (dx0*dx0 + dy0*dy0 + dz0*dz0)^0.5;
if strcmp(mode, 'relative')
    crop = crop*dr;
    dx = crop*dx0/dr;
    dy = crop*dy0/dr;
    dz = crop*dz0/dr;
elseif strcmp(mode, 'absolute')
    crop = crop/dr;
    dx = crop*dx0;
    dy = crop*dy0;
    dz = crop*dz0;
elseif strcmp(mode, 'all')
    if crop > 1
        return;
    end
    ii=1;
	step = crop;
    while (step*ii)<1
        crop = (step*ii)*dr;
        dx(ii) = crop*dx0/dr;
        dy(ii) = crop*dy0/dr;
        dz(ii) = crop*dz0/dr;
        ii=ii+1;
    end
end

for ii=1:length(dx)
    p3(ii,1) = p1(1)-dx(ii);
    p3(ii,2) = p1(2)-dy(ii);
    p3(ii,3) = p1(3)-dz(ii);
end

% ------------------------------------------------------------------------
function [gapRelative] = set_line_crop_AtlasViewerGUI(p1, p2, gridsize)
gapRelative = .05;

if isempty(p1) || isempty(p2)
    return;
end

% Set desired absolute gap distance
gapAbsolute  = 1;
lineSize     = dist3(p1, p2);

% Get relative gap distance
gapRelative = gapAbsolute / lineSize;
% fprintf('Line gap: relative = %0.2f, absolute: %0.2f\n', gapRelative, gapAbsolute);

if gapRelative > .10
    gapRelative = .02;
end

% ---------------------------------------------------------
function [p3] = points_on_line_Homer3(p1, p2, step, mode)

%
% USAGE:
%
%    p3 = points_on_line_Homer3(p1,p2,step,mode)
%
% DESCRIPTION:
%
%    Function takes 2 points in 3d, p1 and p2, and a step length
%    and finds a point on the line formed by extending p1,p2 by the
%    amount (or fraction of p1-p2 distance) step extends away from p1.
%
%    Positive step crops p1,p2 by the step amount.
%    Negative step extends p1,p2 by the step amount.
%
%    mode = 'relative' means step is the distance p3 extends away from p1 relative
%    to the distance between p1 and p2. This is the default if no mode is provided.
%
%    mode = 'absolute' means step is the absolute distance (in unit length) that
%    p3 extends away from p1.
%
%    mode = 'all' returns the set of all points at intervals 1/step on the line
%    segment between p1 and p2.
%
% AUTHOR: Jay Dubb (jdubb@nmr.mgh.harvard.edu)
% DATE:   7/7/2012
%

if ~exist('mode','var')
    mode = 'relative';
end

p3 = p1;

if all(p1==p2)
    return;
end

dx0 = p1(1)-p2(1);
dy0 = p1(2)-p2(2);
dz0 = p1(3)-p2(3);
dr = (dx0*dx0 + dy0*dy0 + dz0*dz0)^0.5;
if strcmp(mode, 'relative')
    crop = step*dr;
    dx = crop*dx0/dr;
    dy = crop*dy0/dr;
    dz = crop*dz0/dr;
elseif strcmp(mode, 'absolute')
    crop = step/dr;
    dx = crop*dx0;
    dy = crop*dy0;
    dz = crop*dz0;
elseif strcmp(mode, 'all')
    if step > 1
        return;
    end
    ii=1;
    while (step*ii)<1
        crop = (step*ii)*dr;
        dx(ii) = crop*dx0/dr;
        dy(ii) = crop*dy0/dr;
        dz(ii) = crop*dz0/dr;
        ii=ii+1;
    end
end

for ii=1:length(dx)
    p3(ii,1) = p1(1)-dx(ii);
    p3(ii,2) = p1(2)-dy(ii);
    p3(ii,3) = p1(3)-dz(ii);
end


text:
%islvector - Test if parameter is a (pseudo) logical vector satisfying an optional list of tests.
%
%  USAGE
%
%    test = islvector(x,test1,test2,...)
%
%    x              parameter to test
%    test1...       optional list of additional tests (see examples below)
%
%  EXAMPLES
%
%    % Test if x is a logical vector
%    islvector(x)
%
%    % Special test: test if x is a logical vector of length 3
%    islvector(x,'#3')
%
%  NOTE
%
%    To be considered logical, the vector should contain only values 0 and 1, but it
%    does not need to actually be of class 'logical' (class(x) could be e.g. 'double').
%
%  SEE ALSO
%
%    See also islscalar, islmatrix, isdmatrix, isdvector, isdscalar, isimatrix, isiscalar,
%    isstring.
%

% Copyright (C) 2010-2011 by Michaël Zugaro
%
% 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.

function test = islvector(x,varargin)

% Check number of parameters
if nargin < 1,
  error('Incorrect number of parameters (type ''help islvector'' for details).');
end

% Test: logical, vector
test = islogical(x) & isvector(x);

% Optional tests
for i = 1:length(varargin),
	try
		if varargin{i}(1) == '#',
			if length(x) ~= str2num(varargin{i}(2:end)), test = false; return; end
		end
	catch err
		error(['Incorrect test ''' varargin{i} ''' (type ''help islvector'' for details).']);
	end
end

function test = islogical(x)

test = builtin('islogical',x) | all(x(:)==0|x(:)==1);


text:
function j = computecost(x, Y, Theta)

n = length(Y); % Number of training examples.

j = 0;

j = (1 / (2 * n)) * sum(((x * Theta) - Y).^2);

end

author: TheAlgorithms

repo: MATLAB-Octave


text:
function H=shadedErrorBar(x,y,errBar,lineProps,transparent)
% function H=shadedErrorBar(x,y,errBar,lineProps,transparent)
%
% Purpose
% Makes a 2-d line plot with a pretty shaded error bar made
% using patch. Error bar color is chosen automatically.
%
% Inputs
% x - vector of x values [optional, can be left empty]
% y - vector of y values or a matrix of n observations by m cases
%     where m has length(x);
% errBar - if a vector we draw symmetric errorbars. If it has a
%          size of [2,length(x)] then we draw asymmetric error bars
%          with row 1 being the upper bar and row 2 being the lower
%          bar. ** alternatively ** errBar can be a cellArray of
%          two function handles. The first defines which statistic
%          the line should be and the second defines the error
%          bar.
% lineProps - [optional,'-k' by default] defines the properties of
%             the data line. e.g.:
%             'or-', or {'-or','markerfacecolor',[1,0.2,0.2]}
% transparent - [optional, 0 by default] if ==1 the shaded error
%               bar is made transparent, which forces the renderer
%               to be openGl. However, if this is saved as .eps the
%               resulting file will contain a raster not a vector
%               image.
%
% Outputs
% H - a structure of handles to the generated plot objects.
%
%
% Examples
% y=randn(30,80); x=1:size(y,2);
% shadedErrorBar(x,mean(y,1),std(y),'g');
% shadedErrorBar(x,y,{@median,@std},{'r-o','markerfacecolor','r'});
% shadedErrorBar([],y,{@median,@std},{'r-o','markerfacecolor','r'});
%
% Overlay two transparent lines
% y=randn(30,80)*10; x=(1:size(y,2))-40;
% shadedErrorBar(x,y,{@mean,@std},'-r',1);
% hold on
% y=ones(30,1)*x; y=y+0.06*y.^2+randn(size(y))*10;
% shadedErrorBar(x,y,{@mean,@std},'-b',1);
% hold off
%
%
% Rob Campbell - November 2009

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Error checking
% error(nargchk(3,5,nargin))

%Process y using function handles if needed to make the error bar
%dynamically
if iscell(errBar) && ~isvector(y)
    fun1=errBar{1};
    fun2=errBar{2};
    errBar=fun2(y);
    y=fun1(y);
elseif ~iscell(errBar) && isvector(y)
    y=y(:)';
else
    error('2nd and 3rd input arguments are not compatible')
end

if isempty(x)
    x=1:length(y);
else
    x=x(:)';
end

if length(x) ~= length(y)
    error('inputs x and y are not of equal lengths')
end

%If only one error bar is specified then we will mirror it, turning it into
%both upper and lower bars.
if length(errBar)==length(errBar(:))
    errBar=repmat(errBar(:)',2,1);
else
    f=find(size(errBar)==2);
    if isempty(f), error('errBar has the wrong size'), end
    if f==2, errBar=errBar'; end
end

if length(x) ~= length(errBar)
    error('inputs x and y must have the same length as errBar')
end

%Set default options
defaultProps={'-k'};
if nargin<4 || isempty(lineProps)
    lineProps=defaultProps;
end
if ~iscell(lineProps)
    lineProps={lineProps};
end

if nargin<5 || ~isnumeric(transparent)
    transparent=0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the main line. We plot this first in order to extract the RGB values
% for the line colour. I am not aware of a function that does this.
H.mainLine=plot(x,y,lineProps{:});

% Work out the color of the shaded region and associated lines
% Using alpha requires the render to be openGL and so you can't
% save a vector image. On the other hand, you need alpha if you're
% overlaying lines. We therefore provide the option of choosing alpha
% or a de-saturated solid colour for the patch surface.

col=get(H.mainLine,'color');
edgeColor=col+(1-col)*0.55;
patchSaturation=0.15; %How de-saturated or transparent to make the patch
if transparent
    faceAlpha=patchSaturation;
    patchColor=col;
    set(gcf,'renderer','openGL')
else
    faceAlpha=1;
    patchColor=col+(1-col)*(1-patchSaturation);
    set(gcf,'renderer','painters')
end

%Calculate the y values at which we will place the error bars
uE=y+errBar(1,:);
lE=y-errBar(2,:);

%Add the error-bar plot elements
holdStatus=ishold;
if ~holdStatus, hold on,  end

%Make the cordinats for the patch
yP=[lE,fliplr(uE)];
xP=[x,fliplr(x)];

%remove any nans otherwise patch won't work
xP(isnan(yP))=[];
yP(isnan(yP))=[];

H.patch=patch(xP,yP,1,'facecolor',patchColor,...
              'edgecolor','none',...
              'facealpha',faceAlpha);

%Make nice edges around the patch.
H.edge(1)=plot(x,lE,'-','color',edgeColor);
H.edge(2)=plot(x,uE,'-','color',edgeColor);

%The main line is now covered by the patch object and was plotted first to
%extract the RGB value of the main plot line. I am not aware of an easy way
%to change the order of plot elements on the graph so we'll just remove it
%and put it back (yuk!)
delete(H.mainLine)
H.mainLine=plot(x,y,lineProps{:});

if ~holdStatus, hold off, end

author: optimaltransport


text:
function M = softMin( D, sigma )
% Calculates the softMin of a vector.
%
% Let D be a vector.  Then the softMin of D is defined as:
%   s = exp(-D/sigma^2) / sum( exp(-D/sigma^2) )
% The softMin is a way of taking a dissimilarity (distance) vector D and
% converting it to a similarity vector s, such that sum(s)==1. If D is an
% NxK array, is is treated as N K-dimensional vectors, and the return is
% likewise an NxK array.  This is useful if D is a distance matrix,
% generated by the likes of pdist2.
%
% Note that as sigma->0, softMin's behavior tends toward that of the
% standard min function.  That is the softMin of a vector D has all zeros
% with a single 1 in the location of the smallest value of D. For example,
% "softMin([.2 .4 .1 .3],eps)" returns "[0 0 1 0]".  As sigma->inf, then
% softMin(D,sigma) tends toward "ones(1,n)/n", where n==length(D).
%
% If D contains the squared euclidean distance between a point y and k
% points xi, then there is a probabilistic interpretation for softMin.  If
% we think of the k points representing equal variant gaussians each with
% mean xi and std sigma, then the softMin returns the relative probability
% of y being generated by each gaussian.
%
% USAGE
%  M = softMin( D, sigma )
%
% INPUTS
%  D       - NxK dissimilarity matrix
%  sigma   - controls 'softness' of softMin
%
% OUTPUTS
%  M       - the softMin (indexes into D)
%
% EXAMPLE - 1
%  C = [0 0; 1 0; 0 1; 1 1]; x=[.7,.3; .1 .2];
%  D = pdist2( x, C ), M = softMin( D, .25 )
%
% EXAMPLE - 2
%  fplot( 'softMin( [0.5 0.2 .4], x )', [0 5] );
%  xlabel('sigma'); ylabel('assignments')
%
% See also PDIST2, SOFTMAX
%
% Piotr's Image&Video Toolbox      Version 2.0
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( sigma==0 ) % special case, make fast
  [~, inds] = min(D,[],2); [n, k] = size(D);
  M = subsToArray( [(1:n)' inds], ones(n,1), [n k] );

else % general case
  M = exp( -D / sigma^2 );
  M(isinf(M))=1e50;
  sumM = sum( M, 2 );
  sumMzero = (sumM==0);
  if( any(sumMzero) )
    [~, inds] = min(D,[],2); [n, k] = size(D);
    Mhard = subsToArray( [(1:n)' inds], ones(n,1), [n k] );
    M( sumMzero, : ) = Mhard( sumMzero, : );
    sumM = sum( M, 2 );
  end
  M = M ./ sumM( :, ones(1,size(M,2)) );
end

author: zouchuhang

repo: LayoutNet


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: MI versus rotation
%
%   - data                 PETCT, Omega=(0,140)x(0,151), level=4:7, m=[128,128]
%   - viewer               viewImage2D
%   - interpolation        splineInter
%   - distance             MI
%   - transformation       rotation2D
%==============================================================================

clear, close all, help(mfilename);

% setup data, interpolation, transformation, regularization
setup2DPETCTData; level = 6; omega = ML{level}.omega; m = ML{level}.m;

viewImage('reset',viewPara{:},'axis','off');
imgModel('reset','imgModel','splineInter','regularizer','moments','theta',1e0);
[T,R] = imgModel('coefficients',ML{level}.T,ML{level}.R,omega,'out',0);
distance('reset','distance','MI');
center = (omega(2:2:end)-omega(1:2:end))'/2;
trafo('reset','trafo','rotation2D','c',center);
fprintf('%20s : %s\n','viewImage',viewImage);
fprintf('%20s : %s\n','imgModel',imgModel);
fprintf('%20s : %s\n','distance',distance);
fprintf('%20s : %s\n','trafo',trafo);

xc = getCellCenteredGrid(omega,m);
Rc = imgModel(R,omega,xc);

% diffImage = @(Tc) viewImage(128+(Tc-Rc)/2,omega,m);
wc = linspace(-pi/2,pi/2,51);
Dc = zeros(size(wc));

% run the loop over all rotations
for j = 1:length(wc),
  yc = trafo(wc(j),xc);                 % compute transformed grid
  Tc = imgModel(T,omega,yc);               % compute transformed image
  [Dc(j),rc] = distance(Tc,Rc,omega,m); % compute distance

  % visualize
  if j == 1, th = [];
    FAIRfigure(1,'figname',mfilename); clf;
    subplot(1,3,1);  viewImage(Rc,omega,m);            th(1) = title('R');
    subplot(1,3,2);  vh = viewImage(Tc,omega,m);       th(2) = title('T(yc)');
    subplot(1,3,3);  ph = plot(wc(1),Dc(1),'r.','markersize',20);
    th(3) = title(sprintf('%s versus rotation',distance));
    axis([wc(1),wc(end),-inf,inf]); hold on;
    axis('auto y')
    set(th,'fontsize',30);
    FAIRpause;
  else
    set(vh,'cdata',reshape(Tc,m)')
    subplot(1,3,3); set(ph,'visible','off');
    plot(wc(1:j),Dc(1:j),'k-','linewidth',2);
    ph = plot(wc(j),Dc(j),'r.','markersize',20);
    drawnow,
    FAIRpause(1/100)
  end;
  fprintf('.'); if ~rem(j,50) || j == length(wc), fprintf('\n'); end;
end;

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


text:
function res = my_non_maxima(I, radii, threshold)

% non maxima supression, odd size mask
mask = 2*radii+1;
max = ordfilt2(I,mask^2,true(mask));

% Maxima and Threshold
res = (I==max) & (I >threshold);

return

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [d dt pred] = bfs(A,u,varargin)
% BFS Compute the breadth first search order.
%
% [d dt pred] = bfs(A,u) returns the distance to each vertex (d) and the
% discover time (dt) in a breadth first search starting from vertex u.
%    d(i) = dt(i) = -1 if vertex i is not reachable from vertex u.
% pred is the predecessor array.  pred(i) = 0 if vertex (i)
% is in a component not reachable from u and i != u.
%
% This method works on directed graphs.
% The runtime is O(V+E).
%
% ... = bfs(A,u,...) takes a set of
% key-value pairs or an options structure.  See set_matlab_bgl_options
% for the standard options.
%   options.target: a special vertex that will stop the search when hit
%       [{'none'} | any vertex number besides the u]
%
% Note: this function does not depend upon the non-zero values of A, but
% only uses the non-zero structure of A.
%
% Example:
%    load graphs/bfs_example.mat
%    d = bfs(A,1)
%
% See also DFS

% David Gleich
% Copyright, Stanford University, 2006-2008

%% History
%  2006-04-19: Initial version
%  2006-05-31: Added full2sparse check
%  2007-04-19: Added target option
%  2008-10-07: Changed options parsing
%%

[trans check full2sparse] = get_matlab_bgl_options(varargin{:});
if full2sparse && ~issparse(A), A = sparse(A); end

options = struct('target', 'none');
options = merge_options(options,varargin{:});

if check, check_matlab_bgl(A,struct()); end

if strcmp(options.target,'none')
    target = 0; % a flag used to denote "no target" to the mex
elseif isa(options.target, 'double')
    target = options.target;
else
    error('matlab_bgl:invalidParameter', ...
        'options.target is not ''none'' or a vertex number.');
end

if (trans)
    A = A';
end

[d dt pred] = bfs_mex(A,u,target);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [prob,sol,fmin] = nls_prob(varargin)
%NLS_PROB  Return an OPTI NLS
%
%   prob = nls_prob(no) return a pre-built optiprob of a saved NLS.
%
%   [prob,sol,fmin] = NLS_prob(no) returns the optimum solution and function
%   eval at the optimum
%
%   no = nls_prob() returns the number of problems available for testing.

%   (C) 2012 Jonathan Currie (I2C2)

% Functions are taken from:
% More, J. J., B. S. Garbow, and K. E. Hillstrom. "Testing Unconstrained
% Optimization Software." ACM Transactions on Mathematical Software 7,
% no. 1 (1981): 17-41.

%Check if just returning no problems
if(nargin < 1)
    prob = 20; sol = []; fmin = [];
    return;
else
    no = varargin{1};
end

%Big switch yard
switch(no)
    case 1
        fun = @(x) [10*(x(2)-x(1)^2); 1 - x(1)];
        ydata = zeros(2,1);
        x0 = [-1.2;1];
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [1,1];
        fmin = 0;

    case 2
        fun = @(x) [-13 + x(1) + ((5 - x(2))*x(2) - 2)*x(2);
                    -29 + x(1) + ((x(2) + 1)*x(2) - 14)*x(2)];
        ydata = zeros(2,1);
        x0 = [0.5;-2];
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [5,4];
        fmin = 0;

    case 3
        fun = @(x) [1e4*x(1)*x(2) - 1;
                    exp(-x(1)) + exp(-x(2)) - 1.0001];
        ydata = zeros(2,1);
        x0 = [0;1];
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [1.098e-5;9.106];
        fmin = 0;

    case 4
        fun = @(x) [x(1) - 1e6; x(2) - 2e-6; x(1)*x(2) - 2];
        ydata = zeros(3,1);
        x0 = [1;1];
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [1e6,2e-6];
        fmin = 0;

    case 5
        fun = @(x) [1.5 - x(1)*(1 - x(2)); 2.25 - x(1)*(1 - x(2)^2); 2.625 - x(1)*(1 - x(2)^3) ];
        ydata = zeros(3,1);
        x0 = [1;1];
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [3;0.5];
        fmin = 0;

    case 6
        i = (1:10)';
        fun = @(x) 2 + 2*i - (exp(0.2578*i) + exp(0.2578*i));
        ydata = zeros(10,1);
        x0 = [0.3;0.4];
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [0.2578;0.2578];
        fmin = 124.36226865;

    case 7
        fun = @more7;
        ydata = zeros(3,1);
        x0 = [-1;0;0];
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [1;0;0];
        fmin = 0;

    case 8
        u = (1:15)';
        v = 16 - u;
        w = zeros(15,1);
        for i = 1:15
            w(i) = min(u(i),v(i));
        end
        y = [0.14;0.18;0.22;0.25;0.29;0.32;0.35;0.39;0.37;0.58;0.73;0.96;1.34;2.10;4.39];
        fun = @(x) (x(1) + u ./ (v*x(2) + w*x(3)));
        x0 = [1;1;1];
        prob = optiprob('fun',fun,'ydata',y,'x0',x0);
        sol = [];
        fmin = 8.21487e-3;

    case 9
        i = (1:15)';
        y = [0.0009;0.0044;0.0175;0.0540;0.1295;0.2420;0.3521;0.3989;0.3521;0.2420;0.1295;0.0540;0.0175;0.0044;0.0009];
        t = (8 - i)./2;
        fun = @(x) x(1)*exp((-x(2)*(t - x(3)).^2)./2);
        x0 = [0.4;1;0];
        prob = optiprob('fun',fun,'ydata',y,'x0',x0);
        sol = [];
        fmin = 1.12793e-8;

    case {10,11}
        i = (1:16)';
        y = [34780;28610;23650;19630;16370;13720;11540;9744;8261;7030;6005;5147;4427;3820;3307;2872];
        t = 45 + 5*i;
        fun = @(x) x(1)*exp(x(2)./(t + x(3)));
        x0 = [0.02;4000;250];
        prob = optiprob('fun',fun,'ydata',y,'x0',x0);
        sol = [];
        fmin = 87.9458;

%     case 11 %not working??
%         m = 100;
%         i = [1:m]';
%         t = i./100;
%         y = 25 + (-50 * log(t)).^(2/3);
%         fun = @(x) exp(-(abs(y*100.*i*x(2)).^x(3))./x(1)) - t;
%         x0 = [5;2.5;0.15];
%         ydata = zeros(m,1);
%         prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
%         sol = [50;25;1.5];
%         fmin = 0;

    case 12
        m = 3;
        i = (1:m)';
        t = 0.1*i;
        fun = @(x) exp(-t*x(1)) - exp(-t*x(2)) - x(3)*exp(-t) - exp(-10*t);
        x0 = [0;10;20];
        ydata = zeros(m,1);
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [1;10;1];
        fmin = 0;

    case 13
        fun = @(x) [x(1) + 10*x(2);
                    sqrt(5)*(x(3) - x(4));
                    (x(3) - 2*x(3))^2;
                    sqrt(10)*(x(1)-x(4))^2];
        x0 = [3;-1;0;1];
        ydata = zeros(4,1);
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [0;0;0;0];
        fmin = 0;

    case 14
        fun = @(x) [10*(x(2)-x(1)^2);
                    1 - x(1);
                    sqrt(90)*(x(4)-x(3)^2);
                    1 - x(3);
                    sqrt(10)*(x(2) + x(4) - 2);
                    10^(-0.5)*(x(2)-x(4))];
        x0 = [3;-1;-3;-1];
        ydata = zeros(6,1);
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [1;1;1;1];
        fmin = 0;

    case 15
        y = [0.1957;0.1947;0.1735;0.16;0.0844;0.0627;0.0456;0.0342;0.0323;0.0235;0.0246];
        u = [4;2;1;0.5;0.25;0.167;0.125;0.1;0.0833;0.0714;0.0625];
        fun = @(x) (x(1)*(u.^2 + u*x(2))) ./ (u.^2 + u*x(3) + x(4));
        x0 = [0.25;0.39;0.415;0.39];
        prob = optiprob('fun',fun,'ydata',-y,'x0',x0);
        sol = [];
        fmin = 3.07505e-4;

    case 16
        m = 20;
        i = (1:m)';
        t = i/5;
        fun = @(x) (x(1) + t*x(2) - exp(t)).^2 + (x(3) + x(4)*sin(t) - cos(t)).^2;
        x0 = [25;5;-5;-1];
        ydata = zeros(m,1);
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [];
        fmin = 85822.2016;

    case 17
        m = 33;
        i = (1:m)';
        t = 10*(i - 1);
        y = [0.844;0.908;0.932;0.936;0.925;0.908;0.881;0.850;0.818;0.784;0.751;0.718;0.685;0.658;0.628;0.603;...
             0.58;0.558;0.538;0.522;0.506;0.490;0.478;0.467;0.457;0.448;0.438;0.431;0.424;0.42;0.414;0.411;0.406];
        fun = @(x) (x(1) + x(2)*exp(-t*x(4)) + x(3)*exp(-t*x(5)));
        x0 = [0.5;1.5;-1;0.01;0.02];
        prob = optiprob('fun',fun,'ydata',-y,'x0',x0);
        sol = [];
        fmin = 5.46489e-5;

    case 18
        m = 13;
        i = (1:m)';
        t = 0.1*i;
        y = exp(-t) - 5*exp(-10*t) + 3*exp(-4*t);
        fun = @(x) x(3)*exp(-t*x(1)) - x(4)*exp(-t*x(2)) + x(6)*exp(-t*x(5));
        x0 = [1;2;1;1;1;1];
        prob = optiprob('fun',fun,'ydata',y,'x0',x0);
        sol = [];
        fmin = 0; %5.65565e-3; seems wrong?

    case 19
        m = 65;
        i = (1:m)';
        t = (i-1)/10;
        y = [1.366;1.191;1.112;1.013;0.991;0.885;0.831;0.847;0.786;0.725;0.746;0.679;0.608;0.655;0.616;0.606;0.602;0.626;0.651;0.724;0.649;0.649;...
            0.694;0.644;0.624;0.661;0.612;0.558;0.533;0.495;0.5;0.423;0.395;0.375;0.372;0.391;0.396;0.405;0.428;0.429;0.523;0.562;0.607;0.653;...
            0.672;0.708;0.633;0.668;0.645;0.632;0.591;0.559;0.597;0.625;0.739;0.710;0.729;0.720;0.636;0.581;0.428;0.292;0.162;0.098;0.054];
        fun = @(x) (x(1)*exp(-t*x(5)) + x(2)*exp(-(t-x(9)).^2*x(6)) + x(3)*exp(-(t-x(10)).^2*x(7)) + x(4)*exp(-(t-x(11)).^2*x(8)));
        x0 = [1.3;0.65;0.65;0.7;0.6;3;5;7;2;4.5;5.5];
        prob = optiprob('fun',fun,'ydata',y,'x0',x0);
        sol = [];
        fmin = 4.01377e-2;

    case 20
        m = 31; n = 12;
        fun = @more20;
        x0 = zeros(n,1);
        ydata = zeros(m,1);
        prob = optiprob('fun',fun,'ydata',ydata,'x0',x0);
        sol = [];
        fmin = 4.72238e-10;

    otherwise
        error('Problem not available or not implemented yet');
end

function j = more7(x)

if(x(1) > 0)
    theta = 1/(2*pi)*atan(x(2)/x(1));
else
    theta = 1/(2*pi)*atan(x(2)/x(1)) + 0.5;
end

j(1,1) = 10*(x(3) - 10*theta);
j(2,1) = 10*(x(1)^2 + x(2)^2)^0.5 - 1;
j(3,1) = x(3);

function f = more20(x)
n = length(x);
j = (2:n)';
j2 = (1:n)';
t = 0.1*(1:29)';
f = zeros(31,1);
for k = 1:29
    f(k) = sum( (j-1).*x(2:end).*t(k).^(j-2) ) - sum( x.*t(k).^(j2-1) ) - 1;
end
f(30) = x(1);
f(31) = x(2) - x(1)^2 - 1;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
filename='Bridge_tetrahedra_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'IPOPT'; 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 f = minus( f, g )
% - MINUS. Minus of two CHEBFUN2V.
%   F - G subtracts the CHEBFUN2V F from G componentwise.
%
%   minus(F, G) is called for the syntax f - g.

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

f = plus( f, uminus( g ) );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun