text:
function write_wobj(OBJ,fullfilename)
% Write objects to a Wavefront OBJ file
%
% write_wobj(OBJ,filename);
%
% OBJ struct containing:
%
% OBJ.vertices : Vertices coordinates
% OBJ.vertices_texture: Texture coordinates
% OBJ.vertices_normal : Normal vectors
% OBJ.vertices_point  : Vertice data used for points and lines
% OBJ.material : Parameters from external .MTL file, will contain parameters like
%           newmtl, Ka, Kd, Ks, illum, Ns, map_Ka, map_Kd, map_Ks,
%           example of an entry from the material object:
%       OBJ.material(i).type = newmtl
%       OBJ.material(i).data = 'vase_tex'
% OBJ.objects  : Cell object with all objects in the OBJ file,
%           example of a mesh object:
%       OBJ.objects(i).type='f'
%       OBJ.objects(i).data.vertices: [n x 3 double]
%       OBJ.objects(i).data.texture:  [n x 3 double]
%       OBJ.objects(i).data.normal:   [n x 3 double]
%
% example reading/writing,
%
%   OBJ=read_wobj('examples\example10.obj');
%   write_wobj(OBJ,'test.obj');
%
% example isosurface to obj-file,
%
%   % Load MRI scan
%   load('mri','D'); D=smooth3(squeeze(D));
%   % Make iso-surface (Mesh) of skin
%   FV=isosurface(D,1);
%   % Calculate Iso-Normals of the surface
%   N=isonormals(D,FV.vertices);
%   L=sqrt(N(:,1).^2+N(:,2).^2+N(:,3).^2)+eps;
%   N(:,1)=N(:,1)./L; N(:,2)=N(:,2)./L; N(:,3)=N(:,3)./L;
%   % Display the iso-surface
%   figure, patch(FV,'facecolor',[1 0 0],'edgecolor','none'); view(3);camlight
%   % Invert Face rotation
%   FV.faces=[FV.faces(:,3) FV.faces(:,2) FV.faces(:,1)];
%
%   % Make a material structure
%   material(1).type='newmtl';
%   material(1).data='skin';
%   material(2).type='Ka';
%   material(2).data=[0.8 0.4 0.4];
%   material(3).type='Kd';
%   material(3).data=[0.8 0.4 0.4];
%   material(4).type='Ks';
%   material(4).data=[1 1 1];
%   material(5).type='illum';
%   material(5).data=2;
%   material(6).type='Ns';
%   material(6).data=27;
%
%   % Make OBJ structure
%   clear OBJ
%   OBJ.vertices = FV.vertices;
%   OBJ.vertices_normal = N;
%   OBJ.material = material;
%   OBJ.objects(1).type='g';
%   OBJ.objects(1).data='skin';
%   OBJ.objects(2).type='usemtl';
%   OBJ.objects(2).data='skin';
%   OBJ.objects(3).type='f';
%   OBJ.objects(3).data.vertices=FV.faces;
%   OBJ.objects(3).data.normal=FV.faces;
%   write_wobj(OBJ,'skinMRI.obj');
%
% Function is written by D.Kroon University of Twente (June 2010)

if(exist('fullfilename','var')==0)
    [filename, filefolder] = uiputfile('*.obj', 'Write obj-file');
    fullfilename = [filefolder filename];
end
[filefolder,filename] = fileparts( fullfilename);

comments=cell(1,4);
comments{1}=' Produced by Matlab Write Wobj exporter ';
comments{2}='';

fid = fopen(fullfilename,'w');
write_comment(fid,comments);

if(isfield(OBJ,'material')&&~isempty(OBJ.material))
    filename_mtl=fullfile(filefolder,[filename '.mtl']);
    fprintf(fid,'mtllib %s\n',filename_mtl);
    write_MTL_file(filename_mtl,OBJ.material)

end

if(isfield(OBJ,'vertices')&&~isempty(OBJ.vertices))
    write_vertices(fid,OBJ.vertices,'v');
end

if(isfield(OBJ,'vertices_point')&&~isempty(OBJ.vertices_point))
    write_vertices(fid,OBJ.vertices_point,'vp');
end

if(isfield(OBJ,'vertices_normal')&&~isempty(OBJ.vertices_normal))
    write_vertices(fid,OBJ.vertices_normal,'vn');
end

if(isfield(OBJ,'vertices_texture')&&~isempty(OBJ.vertices_texture))
    write_vertices(fid,OBJ.vertices_texture,'vt');
end

for i=1:length(OBJ.objects)
    type=OBJ.objects(i).type;
    data=OBJ.objects(i).data;
    switch(type)
        case 'usemtl'
            fprintf(fid,'usemtl %s\n',data);
        case 'f'
            check1=(isfield(OBJ,'vertices_texture')&&~isempty(OBJ.vertices_texture));
            check2=(isfield(OBJ,'vertices_normal')&&~isempty(OBJ.vertices_normal));
            if(check1&&check2)
                for j=1:size(data.vertices,1)
                    fprintf(fid,'f %d/%d/%d',data.vertices(j,1),data.texture(j,1),data.normal(j,1));
                    fprintf(fid,' %d/%d/%d', data.vertices(j,2),data.texture(j,2),data.normal(j,2));
                    fprintf(fid,' %d/%d/%d\n', data.vertices(j,3),data.texture(j,3),data.normal(j,3));
                end
            elseif(check1)
                for j=1:size(data.vertices,1)
                    fprintf(fid,'f %d/%d',data.vertices(j,1),data.texture(j,1));
                    fprintf(fid,' %d/%d', data.vertices(j,2),data.texture(j,2));
                    fprintf(fid,' %d/%d\n', data.vertices(j,3),data.texture(j,3));
                end
            elseif(check2)
                for j=1:size(data.vertices,1)
                    fprintf(fid,'f %d//%d',data.vertices(j,1),data.normal(j,1));
                    fprintf(fid,' %d//%d', data.vertices(j,2),data.normal(j,2));
                    fprintf(fid,' %d//%d\n', data.vertices(j,3),data.normal(j,3));
                end
            else
                for j=1:size(data.vertices,1)
                    fprintf(fid,'f %d %d %d\n',data.vertices(j,1),data.vertices(j,2),data.vertices(j,3));
                end
            end
        otherwise
            fprintf(fid,'%s ',type);
            if(iscell(data))
                for j=1:length(data)
                    if(ischar(data{j}))
                        fprintf(fid,'%s ',data{j});
                    else
                        fprintf(fid,'%0.5g ',data{j});
                    end
                end
            elseif(ischar(data))
                 fprintf(fid,'%s ',data);
            else
                for j=1:length(data)
                    fprintf(fid,'%0.5g ',data(j));
                end
            end
            fprintf(fid,'\n');
    end
end
fclose(fid);

function write_MTL_file(filename,material)
fid = fopen(filename,'w');
comments=cell(1,2);
comments{1}=' Produced by Matlab Write Wobj exporter ';
comments{2}='';
write_comment(fid,comments);

for i=1:length(material)
    type=material(i).type;
    data=material(i).data;
    switch(type)
        case('newmtl')
            fprintf(fid,'%s ',type);
            fprintf(fid,'%s\n',data);
        case{'Ka','Kd','Ks'}
            fprintf(fid,'%s ',type);
            fprintf(fid,'%5.5f %5.5f %5.5f\n',data);
        case('illum')
            fprintf(fid,'%s ',type);
            fprintf(fid,'%d\n',data);
        case {'Ns','Tr','d'}
            fprintf(fid,'%s ',type);
            fprintf(fid,'%5.5f\n',data);
        otherwise
            fprintf(fid,'%s ',type);
            if(iscell(data))
                for j=1:length(data)
                    if(ischar(data{j}))
                        fprintf(fid,'%s ',data{j});
                    else
                        fprintf(fid,'%0.5g ',data{j});
                    end
                end
            elseif(ischar(data))
                fprintf(fid,'%s ',data);
            else
                for j=1:length(data)
                    fprintf(fid,'%0.5g ',data(j));
                end
            end
            fprintf(fid,'\n');
    end
end

comments=cell(1,2);
comments{1}='';
comments{2}=' EOF';
write_comment(fid,comments);
fclose(fid);

function write_comment(fid,comments)
for i=1:length(comments), fprintf(fid,'# %s\n',comments{i}); end

function write_vertices(fid,V,type)
switch size(V,2)
    case 1
        for i=1:size(V,1)
            fprintf(fid,'%s %5.5f\n', type, V(i,1));
        end
    case 2
        for i=1:size(V,1)
            fprintf(fid,'%s %5.5f %5.5f\n', type, V(i,1), V(i,2));
        end
    case 3
        for i=1:size(V,1)
            fprintf(fid,'%s %5.5f %5.5f %5.5f\n', type, V(i,1), V(i,2), V(i,3));
        end
    otherwise
end
switch(type)
    case 'v'
        fprintf(fid,'# %d vertices \n', size(V,1));
    case 'vt'
        fprintf(fid,'# %d texture verticies \n', size(V,1));
    case 'vn'
        fprintf(fid,'# %d normals \n', size(V,1));
    otherwise
        fprintf(fid,'# %d\n', size(V,1));

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function P = spm_imatrix(M)
% returns the parameters for creating an affine transformation
% FORMAT P = spm_imatrix(M)
% M      - Affine transformation matrix
% P      - Parameters (see spm_matrix for definitions)
%___________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner & Stefan Kiebel
% $Id$

% Translations and zooms
%-----------------------------------------------------------------------
R         = M(1:3,1:3);
C         = chol(R'*R);
P         = [M(1:3,4)' 0 0 0  diag(C)'  0 0 0];
if det(R)<0, P(7)=-P(7);end % Fix for -ve determinants

% Shears
%-----------------------------------------------------------------------
C         = diag(diag(C))\C;
P(10:12)  = C([4 7 8]);
R0        = spm_matrix([0 0 0  0 0 0 P(7:12)]);
R0        = R0(1:3,1:3);
R1        = R/R0;

% This just leaves rotations in matrix R1
%-----------------------------------------------------------------------
%[          c5*c6,           c5*s6, s5]
%[-s4*s5*c6-c4*s6, -s4*s5*s6+c4*c6, s4*c5]
%[-c4*s5*c6+s4*s6, -c4*s5*s6-s4*c6, c4*c5]

P(5) = asin(rang(R1(1,3)));
if (abs(P(5))-pi/2)^2 < 1e-9,
    P(4) = 0;
    P(6) = atan2(-rang(R1(2,1)), rang(-R1(3,1)/R1(1,3)));
else
    c    = cos(P(5));
    P(4) = atan2(rang(R1(2,3)/c), rang(R1(3,3)/c));
    P(6) = atan2(rang(R1(1,2)/c), rang(R1(1,1)/c));
end;
return;

% There may be slight rounding errors making b>1 or b<-1.
function a = rang(b)
a = min(max(b, -1), 1);
return;


text:
function value = c1_jac_monomial_integral ( alpha, beta, expon )

%*****************************************************************************80
%
%% C1_JAC_MONOMIAL_INTEGRAL: integral of a monomial with Jacobi weight over C1.
%
%  Discussion:
%
%    value = integral ( -1 <= x <= +1 ) x^expon (1-x)^alpha (1+x)^beta dx
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, the exponent of (1-X) in the weight factor.
%
%    Input, real BETA, the exponent of (1+X) in the weight factor.
%
%    Input, integer EXPON, the exponent.
%
%    Output, real VALUE, the value of the integral.
%
  c = expon;

  if ( mod ( expon, 2 ) == 0 )
    s = +1.0;
  else
    s = -1.0;
  end

  arg1 = - alpha;;
  arg2 =   1.0 + c;
  arg3 =   2.0 + beta + c;
  arg4 = - 1.0;

  value1 = r8_hyper_2f1 ( arg1, arg2, arg3, arg4 );

  arg1 = - beta;
  arg2 =   1.0 + c;
  arg3 =   2.0 + alpha + c;
  arg4 = - 1.0;

  value2 = r8_hyper_2f1 ( arg1, arg2, arg3, arg4 );

  value = gamma ( 1.0 + c ) ...
    * ( s * gamma ( 1.0 + beta  ) * value1 / gamma ( 2.0 + beta  + c ) ...
    +       gamma ( 1.0 + alpha ) * value2 / gamma ( 2.0 + alpha + c ) );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% File c:\M_files\shortcuts\rtdbimod.m
% mca of rtd circuit using pre-screened (bi-modal) inputs;
% uses MATLAB function G2a.m
% Revised and updated 3/10/04
% Since this is a DC circuit, SS arrays A,B,D,E are not used.
clc;clear
% Component values
R1=4.53;R2=34.8;R3=132;R4=9.09;R5=9.09;E1=5;
R6=4.53;R7=27.4;R8=20;R9=20;RT=1.915;
Nom=[R1 R2 R3 R4 R5 R6 R7 R8 R9 RT E1];
Vo=G2a(Nom);
%
% "Real world" tolerances
%
Tinit=0.001;Tlife=0.002;ppm=1e-6;
TC1=50*ppm;TC2=25*ppm;
Thi=Tinit+Tlife+35*TC1;Tlo=-Tinit-Tlife-80*TC1;
Trhi=8.1*1e-4;Trlo=-Trhi;Trefhi=0.02+35*TC2;
Treflo=-0.02-80*TC2;
%
p=1:9;T(1,p)=Tlo;T(2,p)=Thi;
T(1,10)=Trlo;T(2,10)=Trhi;T(1,11)=Treflo;T(2,11)=Trefhi;
%
Nk=10000; % Number of Monte Carlo samples
%
Nc=length(Nom); % Number of components
nb=30; % Number of bins in histograms
Ng=50; % Number of point in ideal Gaussian curve
randn('state',sum(100*clock)); % randomize normal RNG seed
Yn=zeros(Nk,Nc);sp=0.5; % sp is 1/2 of gap width
Rn=zeros(Nk,Nc);
% Get tolerance constants independent of samples k
p=1:Nc;tr1=(T(2,:)-T(1,:))/6;tr2=T(1,:)+1;
%
for w=1:Nc
	k=0;
	while k < Nk
		z=randn;
		if (z<-sp)|(z>sp) % accept only rv's outside of -sp to +sp gap
			k=k+1; % next rv
			Yn(k,w)=z; % store in Yn array
		end
	end
end
for k=1:Nk
   Rn(k,:)=Nom.*(tr1.*(Yn(k,:)+3)+tr2);
   Vm(k)=G2a(Rn(k,:));
end
%
% get Nc input histograms
%
for p=1:Nc
	Vav(p)=mean(Yn(:,p));
	Vsd(p)=3*std(Yn(:,p));
	hin(p,:)=hist(Yn(:,p),nb)/Nk;
	VL1(p)=min(Yn(:,p));VH1(p)=max(Yn(:,p));
	intv1(p)=(VH1(p)-VL1(p))/nb;
	q=1:nb;bin1(p,q)=VL1(p)+intv1(p)*(q-1);
end
%
% get output histogram
%
Vs=std(Vm);Vavg=mean(Vm);
hout=hist(Vm,nb)/Nk;VL2=min(Vm);VH2=max(Vm);
intv2=(VH2-VL2)/nb;
q=1:nb;bin2(q)=VL2+intv2*(q-1);
% Ideal Gaussian curve
intvn2=(VH2-VL2)/Ng;
c1=intv2/(Vs*sqrt(2*pi));
for q=1:Ng
   x1(q)=intvn2*(q-1)+VL2;
   y1(q)=c1*exp((-(x1(q)-Vavg)^2/(2*Vs^2)));
end
%
Vhi2=Vavg+Vs;Vlo2=Vavg-Vs;
Vsr=sprintf('%2.3f\n',3*Vs);Vavgr=sprintf('%2.3f\n',Vavg);
%
subplot(2,1,1)
set(gca,'FontSize',8);
bar(bin1(1,:),hin(1,:),1,'y');
%stairs(bin1(1,:),hin(1,:),'k');
%hold on
%stairs(bin1(2,:),hin(2,:),'b');stairs(bin1(3,:),hin(3,:),'g');
%stairs(bin1(4,:),hin(4,:),'k');stairs(bin1(5,:),hin(5,:),'k');
%stairs(bin1(6,:),hin(6,:),'k');stairs(bin1(7,:),hin(7,:),'k');
%stairs(bin1(8,:),hin(8,:),'k');hold off;
%title('8 of 11 Pre-screened inputs');
title('1 of 11 Pre-screened inputs');
grid off
xlabel('Sigma');
axis([-4 4 0 0.2]);
%
subplot(2,1,2)
bar(bin2,hout,1,'y');
set(gca,'FontSize',[8]);
hold on
h=plot(x1-intv2/2,y1,'k');
hold off
title('RTD output');xlabel('Volts dc')
xlabel('Volts DC');
axis([4.0 4.6 0 0.15]);
text(4.1,0.1,['Vavg=',Vavgr],'FontSize',8);
text(4.1,0.08,['3s=',Vsr],'FontSize',8);
text(4.1,0.06,['Nk = ',num2str(Nk)],'FontSize',8);
%
figure(1)

author: Sable

repo: mcbench-benchmarks


text:
function area = polygonArea3d(poly, varargin)
%POLYGONAREA3D Area of a 3D polygon.
%
%   AREA = polygonArea3d(POLY)
%   POLY is given as a N-by-3 array of vertex coordinates. The resulting
%   area is positive.
%   Works also for polygons given as a cell array of polygons.
%
%   Example
%     % area of a simple 3D square
%     poly = [10 30 20;20 30 20;20 40 20;10 40 20];
%     polygonArea3d(poly)
%     ans =
%        100
%
%     % Area of a 3D mesh
%     [v f] = createCubeOctahedron;
%     polygons = meshFacePolygons(v, f);
%     areas = polygonArea3d(polygons);
%     sum(areas)
%     ans =
%         18.9282
%
%   See also
%     polygons3d, triangleArea3d, polygonArea, polygonCentroid3d

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2012-02-24,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2012 INRA - Cepia Software Platform.

%   HISTORY
%   2013-08-20 add support for multiple polygons

% Check multiple polygons
if iscell(poly) || sum(sum(isnan(poly))) > 0
    % split the polygons into a cell array
    polygons = splitPolygons3d(poly);
    nPolys = length(polygons);

    % compute area of each polygon
    area = zeros(nPolys, 1);
    for  i = 1:nPolys
        area(i) = polygonArea3d(polygons{i});
    end

    return;
end

% put the first vertex at origin (reducing computation errors for polygons
% far from origin)
v0 = poly(1, :);
poly = bsxfun(@minus, poly, v0);

% indices of next vertices
N = size(poly, 1);
iNext = [2:N 1];

% compute cross-product of each elementary triangle formed by origin and
% two consecutive vertices
cp = cross(poly, poly(iNext,:), 2);

% choose one of the triangles as reference for the normal direction
vn = vectorNorm3d(cp);
[tmp, ind] = max(vn); %#ok<ASGLU>
cpRef = cp(ind,:);

% compute the sign of the area of each triangle
% (need to compute the sign explicitely, as the norm of the cross product
% does not keep orientation within supporting plane)
sign_i = sign(dot(cp, repmat(cpRef, N, 1), 2));

% compute area of each triangle, using sign correction
area_i = vectorNorm3d(cp) .* sign_i;

% sum up individual triangles area
area = sum(area_i) / 2;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% Online Bayesian model selection demo.

% We generate data from the model A->B
% and compute the posterior prob of all 3 dags on 2 nodes:
%  (1) A B,  (2) A <- B , (3) A -> B
% Models 2 and 3 are Markov equivalent, and therefore indistinguishable from
% observational data alone.

% We control the dependence of B on A by setting
% P(B|A) = 0.5 - epislon and vary epsilon
% as in Koller & Friedman book p512

% ground truth
N = 2;
dag = zeros(N);
A = 1; B = 2;
dag(A,B) = 1;

ntrials = 100;
ns = 2*ones(1,N);
true_bnet = mk_bnet(dag, ns);
true_bnet.CPD{1} = tabular_CPD(true_bnet, 1, [0.5 0.5]);

% hypothesis space
G = mk_all_dags(N);
nhyp = length(G);
hyp_bnet = cell(1, nhyp);
for h=1:nhyp
  hyp_bnet{h} = mk_bnet(G{h}, ns);
  for i=1:N
    % We must set the CPTs to the mean of the prior for sequential log_marg_lik to be correct
    % The BDeu prior is score equivalent, so models 2,3 will be indistinguishable.
    % The uniform Dirichlet prior is not score equivalent...
    fam = family(G{h}, i);
    hyp_bnet{h}.CPD{i}= tabular_CPD(hyp_bnet{h}, i, 'prior_type', 'dirichlet', ...
				    'CPT', 'unif');
  end
end

clf
seeds = 1:3;
expt = 1;
for seedi=1:length(seeds)
  seed = seeds(seedi);
  rand('state', seed);
  randn('state', seed);

  es = [0.05 0.1 0.15 0.2];
  for ei=1:length(es)
    e = es(ei);
    true_bnet.CPD{2} = tabular_CPD(true_bnet, 2, [0.5+e 0.5-e; 0.5-e 0.5+e]);

    prior = normalise(ones(1, nhyp));
    hyp_w = zeros(ntrials+1, nhyp);
    hyp_w(1,:) = prior(:)';
    LL = zeros(1, nhyp);
    ll = zeros(1, nhyp);
    for t=1:ntrials
      ev = cell2num(sample_bnet(true_bnet));
      for i=1:nhyp
	ll(i) = log_marg_lik_complete(hyp_bnet{i}, ev);
	hyp_bnet{i} = bayes_update_params(hyp_bnet{i}, ev);
      end
      prior = normalise(prior .* exp(ll));
      LL = LL + ll;
      hyp_w(t+1,:) = prior;
    end

    % Plot posterior model probabilities
    % Red = model 1 (no arcs), blue/green = models 2/3 (1 arc)
    % Blue = model 2 (2->1)
    % Green = model 3 (1->2, "ground truth")

    subplot2(length(seeds), length(es), seedi, ei);
    m = size(hyp_w,1);
    h=plot(1:m, hyp_w(:,1), 'r-',  1:m, hyp_w(:,2), 'b-.', 1:m, hyp_w(:,3), 'g:');
    axis([0 m   0 1])
    %title('model posterior vs. time')
    title(sprintf('e=%3.2f, seed=%d', e, seed));
    drawnow
    expt = expt + 1;
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function h=hmatrix3d(p,xx,yy,zz,dd,hh,varargin)

%   Copyright (C) 2004-2006 Per-Olof Persson. See COPYRIGHT.TXT for details.

h=interpn(xx,yy,zz,hh,p(:,1),p(:,2),p(:,3),'*linear');

repo: ifem


text:
function [CV,CE,I,CN,GD,D] = polygonize(V,F,fun,varargin)
  % Polygonize (contour) an implicit function in the spirit of "An Implicit
  % Surface Polygonizer" [Bloomenthal 1994]
  %
  % [CV,CE] = polygonize(V,F,fun)
  %
  % Inputs:
  %   V  #V by dim list of vertex positions of original mesh
  %   F  #F by dim+1 list of element indices into V
  %   fun  function handle so that zero is the desired level set
  % Outputs:
  %   CV  #CV by dim list of contour mesh vertices
  %   CE  #CE by dim list of facet indices into CV
  %   I  #CE list of indices into F
  %   CN  #CN by dim list of unit normal vectors at vertices
  %

  % default values
  delta = [];
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Delta'}, ...
    {'delta'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous funtion to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  flipped_order = flipped_tet_orders();

  max_iters = 10;
  D = fun(V);
  interval = @(DF) any(DF>0.0,2) & any(DF<=0.0,2);
  crosses = @(DE) ...
    (DE(:,1)>0 & DE(:,2)<= 0)*-1 + ...
    (DE(:,2)>0 & DE(:,1)<= 0)*1;
  crossing_F = find(interval(D(F)));
  FF = F(crossing_F,:);
  % simplex size
  ss = size(F,2);
  switch ss
  case 3
    allE = [FF(:,[2 3]);FF(:,[3 1]);FF(:,[1 2])];
  case 4
    allE = ...
      [FF(:,1) FF(:,2); ...
       FF(:,1) FF(:,3); ...
       FF(:,1) FF(:,4); ...
       FF(:,2) FF(:,3); ...
       FF(:,2) FF(:,4); ...
       FF(:,3) FF(:,4) ...
       ];
  end
  sE = sort(allE,2);
  [E,~,EMAP] = unique(sE,'rows');
  cross_dir = crosses(D(E));
  crossing = cross_dir ~= 0;
  J = (1:size(E,1))';
  EE = E(crossing,:);
  J(crossing) = 1:size(EE,1);
  % Blasphemy
  switch ss
  case 3
    CE = sort(reshape(crossing(EMAP),[],3).*reshape(J(EMAP),[],3),2);
    CE = CE(:,2:end);
  case 4
    % CE(f,i) = 0 if ith edge of element f does not cross, otherwise
    % CE(f,i) is the index of the unique edge that does cross
    CE = reshape(crossing(EMAP),[],6).*reshape(J(EMAP),[],6);
    % If 3 edges cross then we can surface with a single triangle
    crossE = reshape(cross_dir(EMAP),[],6);
    IT = find(sum(CE>0,2)==3);
    crossT = crossE(IT,:);
    [CT,ST] = sort(CE(IT,:),2);
    CT = CT(:,end-2:end);
    flippedT = [
      -1 -1 -1 0 0 0 4 5 6 1 3 2
      -1 -1 -1 0 0 0 4 5 6 2 1 3
      -1 0 0 1 1 0 2 3 6 1 5 4
      0 -1 0 -1 0 -1 1 3 5 2 4 6
      0 -1 0 -1 0 1 1 3 5 2 4 6
      0 -1 0 1 0 1 1 3 5 2 4 6
      0 0 -1 0 -1 -1 1 2 4 3 6 5
      0 0 -1 0 1 -1 1 2 4 3 6 5
      0 0 1 0 1 -1 1 2 4 3 5 6
      0 0 1 0 1 1 1 2 4 3 5 6
      0 1 0 -1 0 -1 1 3 5 2 6 4
      1 0 0 -1 -1 0 2 3 6 1 4 5
      1 0 0 1 -1 0 2 3 6 1 4 5
      1 0 0 1 1 0 2 3 6 1 4 5
      1 1 1 0 0 0 4 5 6 1 2 3
      1 1 1 0 0 0 4 5 6 2 3 1
    ];
    fT = ismember([crossT ST],flippedT,'rows');
    CT(fT,:) = fliplr(CT(fT,:));

    IQ = find(sum(CE>0,2)==4);
    crossQ = crossE(IQ,:);
    [CQ,SQ] = sort(CE(IQ,:),2);
    CQ = CQ(:,end-3:end);
    flippedQ = [
      -1 -1 0 0 -1 1 3 4 2 1 6 5
      -1 -1 0 0 1 1 3 4 2 1 6 5
      -1 0 -1 -1 0 -1 2 5 1 3 4 6
      -1 0 -1 1 0 -1 2 5 1 3 4 6
      -1 0 -1 1 0 1 2 5 1 3 4 6
      0 -1 -1 -1 -1 0 1 6 3 2 5 4
      0 1 1 -1 -1 0 1 6 2 3 4 5
      0 1 1 -1 1 0 1 6 2 3 4 5
      0 1 1 1 1 0 1 6 2 3 4 5
      1 0 1 1 0 1 2 5 3 1 6 4
      1 1 0 0 -1 -1 3 4 1 2 5 6
      1 1 0 0 -1 1 3 4 1 2 5 6
    ];
    fQ = ismember([crossQ SQ],flippedQ,'rows');
    %% huh. Need to flip _after_ creating triangles.
    %CQT([fQ;fQ],:) = fliplr(CQT([fQ;fQ],:));
    % or switch the middle two...
    CQ(fQ,:) = CQ(fQ,[1 3 2 4]);
    I = [IT;IQ;IQ];
    I = crossing_F(I);
  end
  % Upper and lower bound on barycenteric coordinate locating =0.5
  EEl= zeros(size(EE,1),1);
  EElV = V(EE(:,1),:);
  Dl = D(EE(:,1));
  EEu= ones(size(EE,1),1);
  Du = D(EE(:,2));
  EEuV = V(EE(:,2),:);
  for iter = 1:max_iters
    EEm = (EEl+EEu)/2;
    CV = 0.5*(EElV+EEuV);
    if iter < max_iters
      Dm = fun(CV);
      front = interval([Dl Dm]);
       EEu(front,:) =  EEm(front,:);
      EEuV(front,:) = CV(front,:);
       Du(front,:) =  Dm(front,:);
       EEl(~front,:) =  EEm(~front,:);
      EElV(~front,:) = CV(~front,:);
       Dl(~front,:) =  Dm(~front,:);
    end
  end

  switch ss
  case 4
    % Flip non-deluanay edges in quads... at least.
    CQT = [CQ(:,[3 1 4]);CQ(:,[2 4 1])];
    A = reshape(internalangles(CV,CQT),[],6);
    % Quads with Non-delaunay diagonals.
    nd = (A(:,1)+A(:,4))>pi;
    CQ(nd,:) = CQ(nd,[3 4 1 2]);
    CQ1 = CQ(:,[3 1 4]);
    CQ2 = CQ(:,[2 4 1]);
    CQ1(nd,:) = fliplr(CQ1(nd,:));
    CQ2(nd,:) = fliplr(CQ2(nd,:));
    CQT = [CQ1;CQ2];
    %% Double check: sum should be zero
    % A = reshape(internalangles(CV,CQT),[],6);
    % % Non-delaunay
    % nd = (A(:,1)+A(:,4))>pi;
    % sum(nd)
    CE = fliplr([CT;CQT]);
    assert(size(CE,2) == ss-1);

    %% How I determined what should be flipped:
    % S = sum(normals(CV,CT).*barycenter(CV,CT),2);
    % A = unique([crossT(S<0,:) ST(S<0,:)],'rows');
    % B = unique([crossT(S>0,:) ST(S>0,:)],'rows');
    % intersect(A,B,'rows');
    % fprintf('    flippedT = [\n');
    % fprintf('      %d %d %d %d %d %d %d %d %d %d %d %d\n',A');
    % fprintf('    ];\n');
    %S = sum(normals(CV,CQ(:,[1 4 3])).*barycenter(CV,CQ(:,[1 4 3])),2);
    %A = unique([crossQ(S<0,:) SQ(S<0,:)],'rows');
    %B = unique([crossQ(S>0,:) SQ(S>0,:)],'rows');
    %intersect(A,B,'rows')
    %fprintf('    flippedQ = [\n');
    %fprintf('      %d %d %d %d %d %d %d %d %d %d %d %d\n',A');
    %fprintf('    ];\n');
  end

  if nargout>=4
    if isempty(delta)
      %  following Bloomenthal's "An Implicit Surface Polygonizer" 1994
      %
      delta = avgedge(V,F)/max_iters^2;
      % More accurate...
      delta = delta*0.0001;
    end
    % Approximate the gradient with finite differences (it'd be cool if fun
    % could also provide the gradient)
    %
    % Combine FD calls into a single call to fun in case fun has a lot of
    % per-call precomputation.
    order = 2
    switch order
    case 1
      % Wasn't computed on the last iteration
      D = fun(CV);
      switch ss
      case 3
        GD = ...
          reshape(fun([CV+[1 0]*delta;CV+[0 1]*delta]),[],2);
      case 4
        GD = reshape( ...
          fun([CV+[1 0 0]*delta;CV+[0 1 0]*delta;CV+[0 0 1]*delta]), ...
          size(CV,1),3);
      end
      % Don't bother dividing by delta (we're normalizing anyway)
      G = GD-D;
    case 2
      delta = delta/2;
      switch ss
      case 3
        GD = ...
          reshape(fun([CV+[1 0]*delta;CV+[0 1]*delta;CV-[1 0]*delta;CV-[0 1]*delta]),[],2,2);
      case 4
        GD = reshape( ...
          fun([ ...
            CV+[1 0 0]*delta;CV+[0 1 0]*delta;CV+[0 0 1]*delta; ...
            CV-[1 0 0]*delta;CV-[0 1 0]*delta;CV-[0 0 1]*delta; ...
            ]), ...
          size(CV,1),3,2);
      end
      G = GD(:,:,1)-GD(:,:,2);
    end

    CN = -normalizerow(G);
  end

end

author: alecjacobson

repo: gptoolbox


text:
function [ t_stop, y_stop ] = p18_stop ( neqn )

%*****************************************************************************80
%
%% P18_STOP returns the stopping point for problem p18.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_STOP, Y_STOP(NEQN), the final data.
%
  y_stop = zeros ( neqn, 1 );

  t_stop = 20.0;
  y_stop = [ ...
   -5.780432953035361e-01; ...
    8.633840009194193e-01; ...
   -9.595083730380727e-01; ...
   -6.504915126712089e-02 ];

  return
end


text:
		% Se sterge spatiul de lucru
clear;
        % Stabilirea numarului de necunoscute
N=1000;
        % Se genereaza matricea coeficientilor (a)
a=rand(N);
		% Se incarca vectorul termenilor liberi b
b=randn(N,1);
        % Setarea cronometrului
t0=cputime;
        % Se rezolva sistemul de ecuatii liniare cu prima metoda prezentata
x1=a\b;
        % Se calculeaza timpul de calcul necesar
t1=cputime-t0;
        % Setarea cronometrului
t0=cputime;
        % Se rezolva sistemul de ecuatii liniare cu a doua metoda prezentata
x2=inv(a)*b;
        % Se calculeaza timpul de calcul necesar
t2=cputime-t0;
        % Se afiseaza timpii de calcul
disp('Timp de calcul metoda 1:'); disp(t1');
disp('Timp de calcul metoda 2:'); disp(t2');


text:
classdef ZXH_CF14 < PROBLEM
% <multi/many> <real> <large/none> <constrained>
% Constrained benchmark MOP proposed by Zhou, Xiang, and He

%------------------------------- Reference --------------------------------
% Y. Zhou, Y. Xiang, and X. He, Constrained multiobjective optimization:
% Test problem construction and performance evaluations, IEEE Transactions
% on Evolutionary Computation, 2021, 25(1): 172-186.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        k;  % Number of constrained variables
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            if isempty(obj.M); obj.M = 2; end
            if isempty(obj.D); obj.D = obj.M+10;  end
            obj.lower    = zeros(1,obj.D) + 1e-10;
            obj.upper    = ones(1,obj.D)  - 1e-10;
            obj.encoding = ones(1,obj.D);
            if obj.M <= 3
                obj.k = obj.M - 1;
            elseif obj.M > 3 && obj.M <= 8
                obj.k = floor(obj.M/2);
            else
                obj.k = 3;
            end
        end
        %% Evaluate multiple solutions
        function Population = Evaluation(obj,varargin)
            PopDec = varargin{1};
            OptX   = 0.2;
            [N,D]  = size(PopDec);
            M      = obj.M;
            % Step 1: Compute cumsum
            Sx = cumsum(PopDec(:,1:M).^2,2,'reverse');
            % Step 2: Compute theta
            THETA = 2/pi*atan(sqrt(Sx(:,2:end))./PopDec(:,1:M-1));
            % Step 3: Calculate Ackley function
            h = 20 - 20 * exp(-0.2 * sqrt(sum((PopDec(:,M+1:end)-OptX).^2,2)/(D-M))) + exp (1) - exp(sum(cos(2 * pi .*(PopDec(:,M+1:end)-OptX)),2)/(D-M));
            % Step 4: Compute T_
            T = (1 - Sx(:,1)).^2 + h;
            % Step 5: Objectives (linear)
            G      = [ones(N,1) cumprod(THETA,2)] .* [1-THETA ones(N,1)];
            PopObj = G .* repmat((1+T),1,M);
            % Step 6: Constraints
            PopCon(:,1) = Sx(:,1) + h - 1;
            PopCon(:,2) = -(Sx(:,1) + h - 1/4);
            for i = 1 : obj.k
                PopCon(:,i+2) = min(min(THETA(:,i)-1/10,4/5-THETA(:,i)),max(2/5-THETA(:,i),THETA(:,i)-7/10));
            end
            Population = SOLUTION(varargin{1},PopObj,PopCon,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = UniformPoint(N,obj.M);
            T = zeros(size(R));
            for i = obj.M-1 : -1 : 1
                T(:,i) = R(:,i+1)./R(:,i)./(1-T(:,i+1)+R(:,i+1)./R(:,i));
            end
            THETA = T(:,1:obj.k);
            Valid = all(THETA<=1/10|THETA>=2/5&THETA<=7/10|THETA>=4/5,2);
            R     = R(Valid,:);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 2
                R = obj.GetOptimum(100);
            elseif obj.M == 3
                a = linspace(0,1,30)';
                R = {a*a',a*(1-a'),(1-a)*ones(size(a'))};
                T2 = R{3}./R{2}./(1+R{3}./R{2});
                T1 = R{2}./(1-T2);
                THETA = cat(3,T1,T2);
                Valid = all(THETA<=1/10|THETA>=2/5&THETA<=7/10|THETA>=4/5,3);
                R{1}(~Valid) = nan;
            else
                R = [];
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [m, a] = polycenter(x,y)
% POLYCENTER  Compute center of mass and area of polygon
%
% [M, A] = POLYCENTER(X, Y)
%
%   M is a 2-vector with the coordinates of the center of mass of a
%   polygon.
%
%   A is the polygon area.
%
%   X, Y are vectors with the coordinates of the polygon vertices.

% Author: Ramon Casero
% Copyright © 2010 University of Oxford
% Version: 0.1.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
error( nargchk( 2, 2, nargin, 'struct' ) );
error( nargoutchk( 0, 2, nargout, 'struct' ) );

% compute polygon area
a = polyarea( x, y );

% close polygon, if open
if (x(1) ~= x(end)) && (y(1) ~= y(end))
    x = x( [1:end 1] );
    y = y( [1:end 1] );
end

% compute x-coordinate of the centroid
m(1) = sum( ...
    ( x( 1:end-1 ) + x( 2:end ) ) .* ( ...
    x( 1:end-1 ) .* y( 2:end ) - x( 2:end ) .* y( 1:end-1 ) ) ...
    ) / 6 / a;

% compute y-coordinate of the centroid
m(2) = sum( ...
    ( y( 1:end-1 ) + y( 2:end ) ) .* ( ...
    x( 1:end-1 ) .* y( 2:end ) - x( 2:end ) .* y( 1:end-1 ) ) ...
    ) / 6 / a;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function x = evaluate_nonlinear(p,x,qq)

% FIX: We have to apply computations to make sure we are evaluating
% expressions such as log(1+sin(x.^2).^2) correctly

if ~isempty(p.bilinears) & all(p.variabletype <= 2) & length(p.evalMap)==0
    x(p.bilinears(:,1)) = x(p.bilinears(:,2)).*x(p.bilinears(:,3));
else
    oldx = 0*p.c;old(1:length(x))=x;
    %try
    x = process_polynomial(x,p);
    %catch
    %    1
    %end
    x = process_evaluations(x,p);
    while norm(x - oldx)>1e-8
        oldx = x;
        x = process_polynomial(x,p);
        x = process_evaluations(x,p);
    end
end

function x = process_evaluations(x,p)
for i = 1:length(p.evalMap)
    arguments = {p.evalMap{i}.fcn,x(p.evalMap{i}.variableIndex)};
    arguments = {arguments{:},p.evalMap{i}.arg{2:end-1}};
    x(p.evalVariables(i)) = feval(arguments{:});
    if ~isempty(p.bilinears)
        x = process_bilinear(x,p);
    end
end

function x = process_bilinear(x,p)
x(p.bilinears(:,1)) = x(p.bilinears(:,2)).*x(p.bilinears(:,3));

function x = process_polynomial(x,p)
x = x(1:length(p.c));
nonlinear = find(p.variabletype);
x(nonlinear) = prod(repmat(x(:)',length(nonlinear),1).^p.monomtable(nonlinear,:),2);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function s = n2s(n)
s = num2str(n);
return

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [grid1, grid2, freq1, freq2, biv, cond2on1, cond1on2, cexp1, cexp2] = bivkernrest(data1,data2,res1,res2,observed)
%[grid1, grid2, freq1, freq2, biv, cond2on1, cond1on2, cexp1, cexp2] = bivkern(data1,data2,res1,res2)
%data1,2 = data for which a kernel density needs to be derived
%res1,2 = step size
%observed = vector of dummies that identifies which elements of data1,2 are
%observations (1) rather than restrictions (0)
%
%restrictions are added by extending the data vectors with
%data1(i) = y for data2(i) = x, observed(i) = 0
%for example, data1(17) = 0, data2(17) = 0, observed(17) = 0
%
%Richard Tol, 13 September 2013

tic

ampl = 0.1; %bandwidth of the restriction
range = 4;  %data are evaluated between minimum minus range*standard deviation and maximum plus range*standard deviation

data1std = std(data1(observed==1));
data1min = res1*round((min(data1(observed==1)) - range*data1std)/res1); %minimum at observed minus range*standard deviation
data1min = max(0,data1min);
data1max = res1*((max(data1(observed==1)) + range*data1std)/res1); %maximum at observed plus range*standard deviation
nopnt = length(data1);
nocon = nopnt - sum(observed);
noobs = nopnt - nocon;

data2std = std(data2(observed==1));
data2min = res2*round((min(data2(observed==1)) - range*data2std)/res2); %minimum at observed minus range*standard deviation
%data2min = max(0,data2min);
data2max = res2*round((max(data2(observed==1)) + range*data2std)/res2); %maximum at observed plus range*standard deviation

%make grid
nogrid1 = round((data1max-data1min)/res1);
grid1 = zeros(nogrid1,1);
grid1(1) = data1min;
for i=2:nogrid1
    grid1(i) = grid1(i-1) + res1;
end

nogrid2 = round((data2max-data2min)/res2);
grid2 = zeros(nogrid2,1);
grid2(1) = data2min;
for i=2:nogrid2
    grid2(i) = grid2(i-1) + res2;
end

%bandwith
sigma = cov(data1(observed==1), data2(observed==1))*1.06^2*noobs^-0.4;
h = inv(sigma);
h1 = 1.06 * noobs^-0.2 * std(data1); %for univariate
h2 = 1.06 * noobs^-0.2 * std(data2); %for univariate

%normal kernel
biv = zeros(nogrid1,nogrid2);
vbiv = biv;
vfreq1 = ones(nogrid1,noobs);
vfreq2 = ones(nogrid2,noobs);
x= zeros(2,1);
for i=1:nopnt,
    %disp(i);
    vfreq1(:,i) = exp(-0.5*((grid1-data1(i))/(h1*observed(i)+ampl*h1*(1-observed(i)))).^2)/(h1*observed(i)+ampl*h1*(1-observed(i)))/sqrt(2*pi);
    vfreq2(:,i) = exp(-0.5*((grid2-data2(i))/(h2*observed(i)+ampl*h2*(1-observed(i)))).^2)/(h2*observed(i)+ampl*h2*(1-observed(i)))/sqrt(2*pi);
    dev1 = grid1-data1(i);
    dev2 = grid2-data2(i);
    for k=1:nogrid1,
        for l=1:nogrid2,
            x1 = dev1(k);
            x2 = dev2(l);
            vbiv(k,l) = exp(-0.5*((h(1,1)*observed(i)+h(1,1)/ampl*(1-observed(i)))*x1^2 + (h(2,2)*observed(i)+h(2,2)/ampl*(1-observed(i)))*x2^2 + (h(1,2)*observed(i)+h(1,2)/ampl*(1-observed(i)))*x1*x2));
        end
    end
    biv = biv + vbiv/sum(sum(vbiv));
end
biv = biv/sum(sum(biv));
freq1 = sum(vfreq1,2);
freq1 = freq1/sum(freq1);
freq2 = sum(vfreq2,2);
freq2 = freq2/sum(freq2);

cond2on1=biv./repmat(sum(biv,2),1,nogrid2);
cond1on2=biv./repmat(sum(biv,1),nogrid1,1);
cond1on2=cond1on2';

cexp1 = sum(cond1on2.*repmat(grid1,1,nogrid2)',2);
cexp2 = sum(cond2on1.*repmat(grid2',nogrid1,1),2);

toc

end


text:
% to prepare the events for inversion based
% on Lu Zhongs code.

report_this_filefun(mfilename('fullpath'));

think
tmp = [a(:,10:12)];
do = ['save ' hodo 'data.inp tmp -ascii'];
err =  ['Error - could not save file ' hodo 'data.inp - permission?'];
err2 = ['errordlg(err);return'];
eval(do,err2)

infi = [hodo 'data.inp'];
outfi = [hodo 'tmpout.dat'];
outfi2 = [hodo 'tmpout2.dat'];

fid = fopen([hodo 'inmifi.dat'],'w');

fprintf(fid,'%s\n',infi);
fprintf(fid,'%s\n',outfi);

fclose(fid);
comm = ['!/bin/rm ' outfi];
eval(comm)

comm = ['!  ' hodi '/stinvers/datasetupSTR_new < ' hodo 'inmifi.dat ' ]
eval(comm)

comm = ['!grep  "1.0" ' outfi  '>'  outfi2];
eval(comm)

comm = ['load ' hodo 'tmpout2.dat'];
eval(comm)

l = a(:,length(a(1,:)));
a(:,10:13) = [];
a = [a  tmpout2  l];

done

author: CelsoReyes


text:
function dataS = populate_planC_structures_field(fieldname, dcmdir_PATIENT_STUDY_SERIES_RTSTRUCTS, scanOriS, attr, ssObj, cObj, optS)
%"populate_planC_structures_field"
%   Given the name of a child field to planC{indexS.structures}, populates that
%   field based on the data contained in dcmdir.PATIENT.STUDY.SERIES.RTSTRUCTS
%   structure passed in, for structure number structNum.
%
%JRA 07/12/06
%YWU Modified 03/01/08
%DK 04/12/09
%   Fixed Coordinate System
%NAV 07/19/16 updated to dcm4che3
%       replaced dcm2ml_Element with getTagValue
%       replaced "get" with "getValue", and
%       "countItems" with "size()"
%
%Usage:
%   dataS = populate_planC_structures_field(fieldname,dcmdir_PATIENT_STUDY_SERIES_RTSTRUCTS);
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

%For easier handling
%global xOffset yOffset;

%STRUCTS = dcmdir_PATIENT_STUDY_SERIES_RTSTRUCTS;

%Default value for undefined fields.
dataS = '';

% if ~exist('attr', 'var')
%     %Grab the dicom object representing the structures.
%     attr = scanfile_mldcm(STRUCTS.file);
% end
%
% %Structure Set ROI Sequence
% %SSRS = attr.getValue(hex2dec('30060020'));
% SSRS = attr.getValue(org.dcm4che3.data.Tag.StructureSetROISequence);
%
% %ROI Contour Sequence
% %RCS = attr.getValue(hex2dec('30060039'));
% RCS = attr.getValue(org.dcm4che3.data.Tag.ROIContourSequence);

% %Count em up.
%
% if ~isempty(RCS)
%     nStructs = RCS.size();
% else
%     nStructs = 0;
% end
% %Structure Set item for this structure.
% ssObj = SSRS.get(structNum - 1);
% %ssObj = org.dcm4che3.data.Attributes(SSRS)
% %ROI Number
% %ROINumber = getTagValue(ssObj, '30060022');
% ROINumber = ssObj.getInts(org.dcm4che3.data.Tag.ROINumber); % vr=IS
%
% %Find the contour object for this structure.
% for s=1:nStructs
%     cObj = RCS.get(s - 1);
%
%     %Referenced ROI Number
%     %RRN = getTagValue(cObj, '30060084');
%     RRN = cObj.getInts(org.dcm4che3.data.Tag.ReferencedROINumber); %IS
%
%     if RRN == ROINumber
%         %We found the correct contour for this structure.
%         break;
%     end
%
% end

switch fieldname

    case 'imageNumber'
        %Currently not implemented

    case 'imageType'
        dataS = 'STRUCTURE';

    case 'caseNumber'
        %Currently not implemented

    case 'roiNumber'
        %dataS = ROINumber;
        dataS = ssObj.getInts(805699618); %org.dcm4che3.data.Tag.ROINumber;

    case 'patientName'
        %Patient's Name
        %nameS = getTagValue(attr, '00100010');
        %dataS = [nameS.FamilyName '^' nameS.GivenName '^' nameS.MiddleName];
        nameObj = javaObject('org.dcm4che3.data.PersonName',attr.getString(1048592)); %org.dcm4che3.data.Tag.PatientName;
        compFamilyName = javaMethod('valueOf','org.dcm4che3.data.PersonName$Component','FamilyName');
        compGivenName = javaMethod('valueOf','org.dcm4che3.data.PersonName$Component','GivenName');
        compMiddleName = javaMethod('valueOf','org.dcm4che3.data.PersonName$Component','MiddleName');
        dataS = [char(nameObj.get(compFamilyName)), '^',...
            char(nameObj.get(compGivenName)), '^',...
            char(nameObj.get(compMiddleName))];

    case 'structureName'
        %ROI Name
        %dataS = getTagValue(ssObj, '30060026');
        dataS = char(ssObj.getString(805699622,0)); %org.dcm4che3.data.Tag.ROIName;

    case 'numberRepresentation'
        %Artifact of RTOG field names, not representative of CERR data.
        dataS = 'CHARACTER';

    case 'structureFormat'
        dataS = 'SCAN-BASED';

    case 'structSetSopInstanceUID'
        %dataS = getTagValue(attr, '00080018');
        dataS = char(attr.getString(524312,0)); %org.dcm4che3.data.Tag.SOPInstanceUID;

    case 'numberOfScans' %aka # of CT slices
        %Referenced Frame of Reference Sequence
        %RFRS = attr.getValue(hex2dec('30060010'));
        %RFRS = attr.getValue(org.dcm4che3.data.Tag.ReferencedFrameOfReferenceSequence);

        %Frame of Reference UID
        %FORUID = getTagValue(ssObj, '30060024');
        %FORUID = char(attr.getStrings(org.dcm4che3.data.Tag.ReferencedFrameOfReferenceUID));

        %Find the series referenced by these contours.  See bottom of file.
        %RSS = getReferencedSeriesSequence(RFRS, FORUID);

        %dataS = getTagValue(RSS,hex2dec('30060016')); % check this
        contourSeq = cObj.getValue(805699648); %org.dcm4che3.data.Tag.ContourSequence;
        if ~isempty(contourSeq)
            dataS = contourSeq.size();
        end

%         %Convert to ML structure format.
%         RSSML = getTagStruct(RSS);
%
%         if ~isempty(RSSML) && ~isempty(RSSML.ContourImageSequence)
%             %# slices in this series.
%             dataS = length(fields(RSSML.ContourImageSequence));
%         else
%             dataS = '';
%         end

    case 'maximumNumberScans'
        %Currently not implemented

    case 'maximumPointsPerSegment'
        %Currently not implemented

    case 'maximumSegmentsPerScan'
        %Currently not implemented

    case 'structureEdition'
        %Currently not implemented

    case 'unitNumber'
        %Currently not implemented

    case 'writer'
        %Manufacturer
        %dataS = getTagValue(attr, '00080070');
        dataS = char(attr.getString(524400,0)); %org.dcm4che3.data.Tag.Manufacturer;

    case 'dateWritten'
        %Structure Set Date.
        %dataS = getTagValue(attr, '30060008');
        dataS = char(attr.getString(805699592,0)); %org.dcm4che3.data.Tag.StructureSetDate;

    case 'structureColor'
        %Currently not implemented

    case 'structureDescription'
        dataS = char(attr.getString(528446,0)); % series description

    case 'roiGenerationAlgorithm'
        dataS = char(attr.getString(805699638,0)); %dataS = getTagValue(ssObj, '30060036');

    case 'roiGenerationDescription'
        dataS = char(attr.getString(805699640,0)); %dataS = getTagValue(ssObj, '30060038');

    case 'studyNumberOfOrigin'
        %Currently not implemented

    case 'contour'
        %Contour Sequence
        if ~cObj.contains(805699648) %org.dcm4che3.data.Tag.ContourSequence %hex2dec('30060040')
            return;
        end
        % cSeq = cObj.getValue(hex2dec('30060040'));
        cSeq = cObj.getValue(805699648); %org.dcm4che3.data.Tag.ContourSequence;

        if cSeq.isEmpty
            return;
        end

        nContours = cSeq.size();

        %optS = opts4Exe([getCERRPath,'CERROptions.json']);
        %contourSliceTol = optS.contourToSliceTolerance;

        for c = 1:nContours
            aContour = cSeq.get(c-1);

            % Referenced SOP instance UID
            % refSeq = aContour.getValue(hex2dec('30060016'));
            refSeq = aContour.getValue(805699606); %org.dcm4che3.data.Tag.ContourImageSequence;
            sopInstanceUID = '';
            sopClassUID = '';
            if ~isempty(refSeq)
                % sopInstanceUID = getTagValue(refSeq.get(0), '00081155');
                sopInstanceUID = char(refSeq.get(0).getString(528725,0)); %org.dcm4che3.data.Tag.ReferencedSOPInstanceUID;
                sopClassUID = char(refSeq.get(0).getString(528720,0)); %SOPClassUID
            end

            dataS(c).sopInstanceUID = sopInstanceUID;
            dataS(c).sopClassUID = sopClassUID;

            %Contour Geometric Type
            %geoType = getTagValue(aContour, '30060042');
            geoType = char(aContour.getString(805699650,0)); %org.dcm4che3.data.Tag.ContourGeometricType;

            switch upper(geoType)
                case 'POINT'
                    warning('Single point contour.');
                case 'OPEN_PLANAR'
                    warning('Open planar contours.')
                case 'OPEN_NONPLANAR'
                    warning('Open, non-planar contours.')
                case 'CLOSED_PLANAR'
                    %Great, continue.
            end

            %Number of Contour Points
            %nPoints = getTagValue(aContour, '30060046');
            nPoints = aContour.getInts(805699654); %org.dcm4che3.data.Tag.NumberOfContourPoints;
            if isempty(nPoints)
                nPoints = 0;
            end
            %             if ~(nPoints > 1)
            %                 dataS(i).segments = [];
            %                 continue;
            %             end

            %Contour Data
            try
                %data    = getTagValue(aContour, '30060050');
                data    = aContour.getDoubles(805699664); %org.dcm4che3.data.Tag.ContourData;
            catch
                disp('vacant contour found ...');
                dataS(c).segments = [];
                continue;
            end
            %Reshape Contour Data
            data    = reshape(data, [3 nPoints])';

            %data(:,3) = -data(:,3); %Z is always negative to match RTOG spec

            % Flip based on pt position

            %Get assoc scanUID
            assocUID = getAssocScanUID(attr,ssObj);
            assocScanNum = strcmp(assocUID,{scanOriS.scanUID});
            if ~any(assocScanNum) && length(scanOriS)==1
                assocScanNum = 1;
            end
            if any(assocScanNum)
                imgOriV = scanOriS(assocScanNum).imageOrientationPatient;
            else
                imgOriV = [1,0,0,0,1,0]'; % HFS
            end
%              if ~isempty(imgOri)
                  data = convertCoordinates(data, imgOriV);
%              else
%                     %Default it to HFS
%                 data(:,3) = -data(:,3); %Z is always negative to match RTOG spec
%                 data(:,2) = -data(:,2);
%              end

            %Convert from DICOM mm to CERR cm.
            data    = data / 10;

            %Replicate the last data point: CERR needs first/last identical
            if ~isempty(data)
                data(end+1,:) = data(1,:);
%% REMOVED IN LATEST CERR
%{
                if length(unique(data(:,3))) > 1;

                    a = abs(diff(unique(data(:,3))));

                    if (max(a) > contourSliceTol)
                        %ROI Name
                        name = getTagValue(ssObj, '30060026');

                        warning(['CERR does not support out-of-plane contours. Skipping contour ' num2str(i) ' in structure ' name '.']);
                        continue;
                    end
                end
%}
            end

            dataS(c).segments = data;

        end

        %Segments are inserted in order of storage in the DICOM sequence,
        %meaning segments must be sorted and rearranged to match the slice
        %order of the associated scan.   This is handled by a function that
        %has access to both the scan and contour data, dcmdir2planC.

    case 'rasterSegments'
        %Implementation not necessary

    case 'DSHPoints'
        %Implementation not necessary

    case 'orientationOfStructure'
        %In CERR, all structures are oriented transversely.
        dataS = 'TRANSVERSE';

    case 'transferProtocol'
        dataS = 'DICOM';

    case 'referencedFrameOfReferenceUID'
        %dataS  = getTagValue(attr, '00200052');
        %dataS = char(attr.getStrings(org.dcm4che3.data.Tag.FrameOfReferenceUID));
        %dataS = char(attr.getString(2097234,0));
        dataS = char(ssObj.getString(805699620,0)); %org.dcm4che3.data.Tag.ReferencedFrameOfReferenceUID;

    case 'referencedSeriesUID'

        dataS = '';

        refFrameOfReferenceSeq = attr.getValue(805699600);
        if ~refFrameOfReferenceSeq.isEmpty
            refFrameOfReferenceSeq = refFrameOfReferenceSeq.get(0);
        else
            return
        end

        reRefStudySequence = refFrameOfReferenceSeq.getValue(805699602);
        if ~reRefStudySequence.isEmpty
            reRefStudySequence = reRefStudySequence.get(0);
        else
            return;
        end

        rtRefSeriesSequence = reRefStudySequence.getValue(805699604);
        if ~rtRefSeriesSequence.isEmpty
            rtRefSeriesSequence = rtRefSeriesSequence.get(0);
        else
            return
        end

        dataS = char(rtRefSeriesSequence.getString(2097166,0));

    case 'DICOMHeaders'
        %Currently not implemented
        %Read all the dcm data into a MATLAB struct.
        if strcmpi(optS.saveDICOMheaderInPlanC,'yes')
            dataS = getTagStruct(attr);
        end
%
%         %Remove contoursequence data to avoid storing huge amounts of redundant data.
%         try
%             dataS = rmfield(dataS, 'ROIContourSequence');
%         end

%         dataS = '';

    case 'visible'
        dataS = true; % default to visible

    case 'associatedScan'
        %Currently not implemented

    case 'strUID'
        dataS = createUID('structure');

    case 'assocScanUID'
        %wy, use the Referenced frame of reference UID to associate structures to scan.
        %SSRS = attr.get(org.dcm4che2.data.Tag.StructureSetROISequence);
        %SSRS_1 = SSRS.getDicomObject(0);
        %dataS = char(SSRS_1.getString(org.dcm4che2.data.Tag.ReferencedFrameofReferenceUID));

        %dataS = getTagValue(SSRS.get(hex2dec('00200052')));
        %dataS = ['CT.',dataS];

        dataS = getAssocScanUID(attr,ssObj);

%         %Convert to ML structure format.
%         RSSML = getTagStruct(RSS);
%
%         if ~isempty(RSSML)
%             %UID of series structures were contoured on.
%             dataS = RSSML.SeriesInstanceUID;
%             dataS = ['CT.',dataS];
%         else
%             dataS = '';
%         end

    otherwise
        %         warning(['DICOM Import has no methods defined for import into the planC{indexS.structures}.' fieldname ' field, leaving empty.']);
end

function attr = getReferencedSeriesSequence(Referenced_Frame_Of_Reference_Sequence, Frame_of_Reference_UID);
%"getReferencedSeriesSequence"
%   Searches the RFRS for a match to the FoRUID of a structure and returns
%   the subfield of the RFRS containing information on the series that the
%   contours were defined on.  Not a general function since it is really
%   only used when populating structure fields.

%Easier reading.
referencedSeq = Referenced_Frame_Of_Reference_Sequence;
refUID = Frame_of_Reference_UID;

if ~isempty(referencedSeq)
    nRFRS = referencedSeq.size();
else
    nRFRS = 0;
end
%Search the RFR sequence for the UID matching the FORUID.
for i=1:nRFRS

    myRFR = referencedSeq.get(i-1);

    % RFRUID = getTagValue(myRFR, '00200052');
    RFRUID = char(myRFR.getString(2097234,0)); %org.dcm4che3.data.Tag.FrameOfReferenceUID;

    if isequal(RFRUID,refUID)

        %RT Referenced Study Sequence
        %RTRSS = myRFR.getValue(hex2dec('30060012'));
        RTRSS = myRFR.getValue(805699602); %org.dcm4che3.data.Tag.RTReferencedStudySequence;

        if isempty(RTRSS)
            break;
        end

        RTRSS_1 = RTRSS.get(0);

        %RT Referenced Series Sequence
        %RTRSS = RTRSS_1.getValue(hex2dec('30060014'));
        RTRSS = RTRSS_1.getValue(805699604); %org.dcm4che3.data.Tag.RTReferencedSeriesSequence;

        attr = RTRSS.get(0);

    end

end

if ~exist('attr', 'var')
    warning('No explicit association between structures and a scan could be found.');
    attr = javaObject('org.dcm4che3.data.Attributes');
end

end %End of function

function assocScanUID = getAssocScanUID(attr,ssObj)
%commented by wy
        %Referenced Frame of Reference Sequence
        %referencedSeq = attr.getValue(hex2dec('30060010'));
        referencedSeq = attr.getValue(805699600); %org.dcm4che3.data.Tag.ReferencedFrameOfReferenceSequence;

        %Frame of Reference UID
        %refUID = getTagValue(ssObj, '30060024');
        refUID = char(ssObj.getString(805699620,0)); %org.dcm4che3.data.Tag.ReferencedFrameOfReferenceUID;

        %Find the series referenced by these contours.  See bottom of file.
        refSerSeq = getReferencedSeriesSequence(referencedSeq, refUID);

        %assocScanUID = getTagValue(refSerSeq,'0020000E');
        assocScanUID = char(refSerSeq.getString(2097166,0)); %org.dcm4che3.data.Tag.SeriesInstanceUID;

        assocScanUID = ['CT.',assocScanUID];

end %End of function

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function data_dir = GetCurrentDataDir()
    data_dir = 'C:\Users\Xiu\Dropbox (Personal)\Janelia2015\GUI_data';
end

author: xiuyechen

repo: FishExplorer


text:
function K = from_data(kern,dat1,dat2,ind1,ind2,kerParam),

    K=get_x(dat1,ind1,dat2.index(ind2))';

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%%%test for Random MultiGraphs
clear all;
clc;

addpath('./rmg/');
addpath('./util/');

%===== basic parameters =================================
bandNum = 4; % number of band for LPE band selection
w = 7;       % patch size
% LBP feature extraction
r = 1;  nr = 8;
% number of graphs
% for computationa efficiency, number of graphs are set as 4.
% you can set more graphs to obtain better performance
kg = 4;

fprintf('... ... loading data begin ...\n');
load IndianPines_Data.mat;
fprintf('... ... loading data finished !!! \n');

% training number for Inidan_Pines dataset
CTrain = [5 143 83 24 48 73 3 48 2 97 246 59 21 127 39 9];

%===== basic parameters ===============================
no_class = max(gth(:));

% data normalization
fprintf(' ... ... data normalization    ... ...\n');
Data = z./max(z(:));
[ylen, xlen, spec_dim] = size(Data);

% band selection
fprintf(' ... ... band selection        ... ...\n');
X = reshape(Data, ylen*xlen, spec_dim);
Psi = PCA_Train(X', bandNum);
X = X*Psi;
DataTmp = reshape(X, ylen, xlen, size(Psi,2));
clear X Psi;

mapping = getmapping(nr,'u2');
fprintf(' ... ... LBP feature extraction begin ... ...\n');
Feature_P = LBP_feature_global(DataTmp, r, nr, mapping, w, gth);
clear nr r z DataTmp;

lbp_dim = size(Feature_P, 3);
% spatial  data
DataSpat = NewScale(reshape(Feature_P, ylen*xlen, lbp_dim));
% spectral data
DataSpec = NewScale(reshape(Data, ylen*xlen, spec_dim));
% spatial and spectral data combination
DataSpec = DataSpec(:, 1:150);
Data_spec_spat = [DataSpat, DataSpec];
clear DataSpat DataSpec Data Feature_P bandNum;
clear lbp_dim mapping w;

Data = []; Labels = [];

for i = 1: no_class
    pos = find(gth==i);
    Data = [Data; Data_spec_spat(pos, :)];
    Labels = [Labels, length(pos)];
end
clear  Data_spec_spat;

DataTrn = []; DataTst = [];  CTest = [];
k = 0;
for i = 1: no_class
    Data_tmp = Data((k+1):(Labels(i)+k), :);
    k = Labels(i) + k;
    index_i = randperm(Labels(i));
    DataTrn = [DataTrn; Data_tmp(index_i(1:CTrain(i)), :)];
    index_i = find(gth==i);
    CTest(i) = length(index_i);
end

DataTst = Data;

clear k Data_tmp Data index_i;

TrnLab = []; TstLab = [];
for jj = 1: length(CTrain)
   TrnLab = [TrnLab; jj * ones(CTrain(jj),1)];
end
for jj = 1: length(CTest)
   TstLab = [TstLab; jj * ones(CTest(jj),1)];
end

% Scale the data
X=[DataTrn;DataTst];
[N,Dim] = size(X);
clear CTest CTrain gth DataTrn DataTst;

fprintf('... ... Graph number:%d ... ...\n', kg);

%kf=floor(log2(Dim)+1);
% gaofeng revised code 2017/04/08
kf = floor(Dim/4);

label_index = find(TrnLab~=0);
labels = [TrnLab;TstLab];
[G,F]  = MultiGraphs(X,labels,label_index,kg,kf);

[val, predict_res]=max(F,[],2);

[Pr, ConfMat] = GetAccuracy(predict_res(length(label_index)+1:end), ...
                            labels(length(label_index)+1:end));

fprintf(' ... ... Final Accuracy: %f\n', Pr.OA);
fprintf(' ... ... Final Kappa: %f\n', Pr.Kappa);


text:
function [p,ellipse]=phantom3d(varargin)

%PHANTOM3D Three-dimensional analogue of MATLAB Shepp-Logan phantom
%   P = PHANTOM3D(DEF,N) generates a 3D head phantom that can
%   be used to test 3-D reconstruction algorithms.
%
%   DEF is a string that specifies the type of head phantom to generate.
%   Valid values are:
%
%      'Shepp-Logan'            A test image used widely by researchers in
%                               tomography
%      'Modified Shepp-Logan'   (default) A variant of the Shepp-Logan phantom
%                               in which the contrast is improved for better
%                               visual perception.
%
%   N is a scalar that specifies the grid size of P.
%   If you omit the argument, N defaults to 64.
%
%   P = PHANTOM3D(E,N) generates a user-defined phantom, where each row
%   of the matrix E specifies an ellipsoid in the image.  E has ten columns,
%   with each column containing a different parameter for the ellipsoids:
%
%     Column 1:  A      the additive intensity value of the ellipsoid
%     Column 2:  a      the length of the x semi-axis of the ellipsoid
%     Column 3:  b      the length of the y semi-axis of the ellipsoid
%     Column 4:  c      the length of the z semi-axis of the ellipsoid
%     Column 5:  x0     the x-coordinate of the center of the ellipsoid
%     Column 6:  y0     the y-coordinate of the center of the ellipsoid
%     Column 7:  z0     the z-coordinate of the center of the ellipsoid
%     Column 8:  phi    phi Euler angle (in degrees) (rotation about z-axis)
%     Column 9:  theta  theta Euler angle (in degrees) (rotation about x-axis)
%     Column 10: psi    psi Euler angle (in degrees) (rotation about z-axis)
%
%   For purposes of generating the phantom, the domains for the x-, y-, and
%   z-axes span [-1,1].  Columns 2 through 7 must be specified in terms
%   of this range.
%
%   [P,E] = PHANTOM3D(...) returns the matrix E used to generate the phantom.
%
%   Class Support
%   -------------
%   All inputs must be of class double.  All outputs are of class double.
%
%   Remarks
%   -------
%   For any given voxel in the output image, the voxel's value is equal to the
%   sum of the additive intensity values of all ellipsoids that the voxel is a
%   part of.  If a voxel is not part of any ellipsoid, its value is 0.
%
%   The additive intensity value A for an ellipsoid can be positive or negative;
%   if it is negative, the ellipsoid will be darker than the surrounding pixels.
%   Note that, depending on the values of A, some voxels may have values outside
%   the range [0,1].
%
%   Example
%   -------
%        ph = phantom3d(128);
%        figure, imshow(squeeze(ph(64,:,:)))
%
%   Copyright 2005 Matthias Christian Schabel (matthias @ stanfordalumni . org)
%   University of Utah Department of Radiology
%   Utah Center for Advanced Imaging Research
%   729 Arapeen Drive
%   Salt Lake City, UT 84108-1218
%
%   This code is released under the Gnu Public License (GPL). For more information,
%   see : http://www.gnu.org/copyleft/gpl.html
%
%   Portions of this code are based on phantom.m, copyrighted by the Mathworks
%

[ellipse,n] = parse_inputs(varargin{:});

p = zeros([n n n]);

rng =  ( (0:n-1)-(n-1)/2 ) / ((n-1)/2);

[x,y,z] = meshgrid(rng,rng,rng);

coord = [flatten(x); flatten(y); flatten(z)];

p = flatten(p);

for k = 1:size(ellipse,1)
   A = ellipse(k,1);            % Amplitude change for this ellipsoid
   asq = ellipse(k,2)^2;        % a^2
   bsq = ellipse(k,3)^2;        % b^2
   csq = ellipse(k,4)^2;        % c^2
   x0 = ellipse(k,5);           % x offset
   y0 = ellipse(k,6);           % y offset
   z0 = ellipse(k,7);           % z offset
   phi = ellipse(k,8)*pi/180;   % first Euler angle in radians
   theta = ellipse(k,9)*pi/180; % second Euler angle in radians
   psi = ellipse(k,10)*pi/180;  % third Euler angle in radians

   cphi = cos(phi);
   sphi = sin(phi);
   ctheta = cos(theta);
   stheta = sin(theta);
   cpsi = cos(psi);
   spsi = sin(psi);

   % Euler rotation matrix
   alpha = [cpsi*cphi-ctheta*sphi*spsi   cpsi*sphi+ctheta*cphi*spsi  spsi*stheta;
            -spsi*cphi-ctheta*sphi*cpsi  -spsi*sphi+ctheta*cphi*cpsi cpsi*stheta;
            stheta*sphi                  -stheta*cphi                ctheta];

   % rotated ellipsoid coordinates
   coordp = alpha*coord;

   idx = find((coordp(1,:)-x0).^2./asq + (coordp(2,:)-y0).^2./bsq + (coordp(3,:)-z0).^2./csq <= 1);
   p(idx) = p(idx) + A;
end

p = reshape(p,[n n n]);

return;

function out = flatten(in)

out = reshape(in,[1 prod(size(in))]);

return;

function [e,n] = parse_inputs(varargin)
%  e is the m-by-10 array which defines ellipsoids
%  n is the size of the phantom brain image

n = 128;     % The default size
e = [];
defaults = {'shepp-logan', 'modified shepp-logan', 'yu-ye-wang'};

for i=1:nargin
   if ischar(varargin{i})         % Look for a default phantom
      def = lower(varargin{i});
      idx = strmatch(def, defaults);
      if isempty(idx)
         eid = sprintf('Images:%s:unknownPhantom',mfilename);
         msg = 'Unknown default phantom selected.';
         error(eid,'%s',msg);
      end
      switch defaults{idx}
      case 'shepp-logan'
         e = shepp_logan;
      case 'modified shepp-logan'
         e = modified_shepp_logan;
      case 'yu-ye-wang'
         e = yu_ye_wang;
      end
   elseif numel(varargin{i})==1
      n = varargin{i};            % a scalar is the image size
   elseif ndims(varargin{i})==2 && size(varargin{i},2)==10
      e = varargin{i};            % user specified phantom
   else
      eid = sprintf('Images:%s:invalidInputArgs',mfilename);
      msg = 'Invalid input arguments.';
      error(eid,'%s',msg);
   end
end

% ellipse is not yet defined
if isempty(e)
   e = modified_shepp_logan;
end

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Default head phantoms:   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function e = shepp_logan

e = modified_shepp_logan;
e(:,1) = [1 -.98 -.02 -.02 .01 .01 .01 .01 .01 .01];

return;

function e = modified_shepp_logan
%
%   This head phantom is the same as the Shepp-Logan except
%   the intensities are changed to yield higher contrast in
%   the image.  Taken from Toft, 199-200.
%
%         A      a     b     c     x0      y0      z0    phi  theta    psi
%        -----------------------------------------------------------------
e =    [  1  .6900  .920  .810      0       0       0      0      0      0
        -.8  .6624  .874  .780      0  -.0184       0      0      0      0
        -.2  .1100  .310  .220    .22       0       0    -18      0     10
        -.2  .1600  .410  .280   -.22       0       0     18      0     10
         .1  .2100  .250  .410      0     .35    -.15      0      0      0
         .1  .0460  .046  .050      0      .1     .25      0      0      0
         .1  .0460  .046  .050      0     -.1     .25      0      0      0
         .1  .0460  .023  .050   -.08   -.605       0      0      0      0
         .1  .0230  .023  .020      0   -.606       0      0      0      0
         .1  .0230  .046  .020    .06   -.605       0      0      0      0 ];

return;

function e = yu_ye_wang
%
%   Yu H, Ye Y, Wang G, Katsevich-Type Algorithms for Variable Radius Spiral Cone-Beam CT
%
%         A      a     b     c     x0      y0      z0    phi  theta    psi
%        -----------------------------------------------------------------
e =    [  1  .6900  .920  .900      0       0       0      0      0      0
        -.8  .6624  .874  .880      0       0       0      0      0      0
        -.2  .4100  .160  .210   -.22       0    -.25    108      0      0
        -.2  .3100  .110  .220    .22       0    -.25     72      0      0
         .2  .2100  .250  .500      0     .35    -.25      0      0      0
         .2  .0460  .046  .046      0      .1    -.25      0      0      0
         .1  .0460  .023  .020   -.08    -.65    -.25      0      0      0
         .1  .0460  .023  .020    .06    -.65    -.25     90      0      0
         .2  .0560  .040  .100    .06   -.105    .625     90      0      0
        -.2  .0560  .056  .100      0    .100    .625      0      0      0 ];

return;


text:
function b = mv_st ( m, n, nz_num, row, col, a, x )

%*****************************************************************************80
%
%% MV_ST multiplies a sparse triple matrix times a vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, integer NZ_NUM, the number of nonzero values.
%
%    Input, integer ROW(NZ_NUM), COL(NZ_NUM), the row and column indices.
%
%    Input, real A(NZ_NUM), the nonzero values in the M by N matrix.
%
%    Input, real X(N), the vector to be multiplied.
%
%    Output, real B(M), the product A*X.
%
  b = zeros ( m, 1 );

  for k = 1 : nz_num
    b(row(k)) = b(row(k)) + a(k) * x(col(k));
  end

  return
end


%==============================================================================
% 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/
%==============================================================================
%
% function [Dc,rc,dD,dr,d2psi] = SSDmex(Tc,Rc,omega,m,varargin)
%
% Wrapper for C file for Sum of Squared Differences based distance measure,
% see also SSD for details
% see also distances/contents
%==============================================================================

function [Dc,rc,dD,dr,d2psi] = SSDmex(Tc,Rc,omega,m,varargin)

if nargin == 0
    help(mfilename);
    setup2DhandData
    xc = getCellCenteredGrid(omega,m);
    Tc = linearInter(dataT,omega,xc);
    Rc = linearInter(dataR,omega,xc);
    D0 = feval(mfilename,Tc,Rc,omega,m);
    fprintf('%s: distance = %s\n',mfilename,num2str(D0))
    return
end

dD = []; dr = []; d2psi = [];
doDerivative = (nargout > 3);

for k=1:2:length(varargin)      % overwrite default parameter
    eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end

hd = prod((omega(2:2:end)-omega(1:2:end))./m); % voxel size for integration

if ~doDerivative
    try
        [Dc,rc] = SSDmexC(Tc,Rc,hd);
    catch err
        FAIRerror(err);
    end
else
    try
        [Dc,rc,dD,dr,d2psi] = SSDmexC(Tc,Rc,hd);
    catch err
        FAIRerror(err);
    end
end

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


text:
%-----------------------------------------------------------------------
% Job saved on 19-Jan-2014 14:08:03 by cfg_util (rev $Rev: 7703 $)
% spm SPM - SPM12b (beta)
% cfg_basicio BasicIO - Unknown
%-----------------------------------------------------------------------
matlabbatch{1}.spm.meeg.source.headmodel.val = 1;
matlabbatch{1}.spm.meeg.source.headmodel.comment = '';
matlabbatch{1}.spm.meeg.source.headmodel.meshing.meshes.template = 1;
matlabbatch{1}.spm.meeg.source.headmodel.meshing.meshres = 2;
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.useheadshape = 0;
matlabbatch{1}.spm.meeg.source.headmodel.forward.eeg = 'EEG BEM';
matlabbatch{1}.spm.meeg.source.headmodel.forward.meg = 'Single Shell';
matlabbatch{2}.spm.tools.beamforming.data.D(1) = cfg_dep('M/EEG head model specification: M/EEG dataset(s) with a forward model', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','D'));
matlabbatch{2}.spm.tools.beamforming.data.val = 1;
matlabbatch{2}.spm.tools.beamforming.data.gradsource = 'inv';
matlabbatch{2}.spm.tools.beamforming.data.space = 'MNI-aligned';
matlabbatch{2}.spm.tools.beamforming.data.overwrite = 1;
matlabbatch{3}.spm.tools.beamforming.sources.BF(1) = cfg_dep('Prepare data: BF.mat file', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{3}.spm.tools.beamforming.sources.reduce_rank = [2 3];
matlabbatch{3}.spm.tools.beamforming.sources.keep3d = 1;
matlabbatch{3}.spm.tools.beamforming.sources.plugin.mesh.orient = 'original';
matlabbatch{3}.spm.tools.beamforming.sources.plugin.mesh.fdownsample = 1;
matlabbatch{3}.spm.tools.beamforming.sources.plugin.mesh.flip = false;
matlabbatch{3}.spm.tools.beamforming.sources.visualise = 1;
matlabbatch{4}.spm.tools.beamforming.features.BF(1) = cfg_dep('Define sources: BF.mat file', substruct('.','val', '{}',{3}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{4}.spm.tools.beamforming.features.whatconditions.all = 1;
matlabbatch{4}.spm.tools.beamforming.features.woi = [-Inf Inf];
matlabbatch{4}.spm.tools.beamforming.features.modality = {'MEG'};
matlabbatch{4}.spm.tools.beamforming.features.fuse = 'no';
matlabbatch{4}.spm.tools.beamforming.features.plugin.cov.foi = [0 Inf];
matlabbatch{4}.spm.tools.beamforming.features.plugin.cov.taper = 'none';
matlabbatch{4}.spm.tools.beamforming.features.regularisation.manual.lambda = 1;
matlabbatch{4}.spm.tools.beamforming.features.bootstrap = false;
matlabbatch{5}.spm.tools.beamforming.inverse.BF(1) = cfg_dep('Covariance features: BF.mat file', substruct('.','val', '{}',{4}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{5}.spm.tools.beamforming.inverse.plugin.eloreta.regularisation = 0.05;
matlabbatch{6}.spm.tools.beamforming.output.BF(1) = cfg_dep('Inverse solution: BF.mat file', substruct('.','val', '{}',{5}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_power.whatconditions.all = 1;
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_power.sametrials = false;
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_power.contrast = 1;
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_power.result = 'singleimage';
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_power.scale = 2;
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_power.powermethod = 'trace';
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_power.modality = 'MEG';
matlabbatch{7}.spm.tools.beamforming.write.BF(1) = cfg_dep('Output: BF.mat file', substruct('.','val', '{}',{6}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{7}.spm.tools.beamforming.write.plugin.gifti.normalise = 'no';
matlabbatch{7}.spm.tools.beamforming.write.plugin.gifti.space = 'mni';
matlabbatch{7}.spm.tools.beamforming.write.plugin.gifti.visualise = 2;

author: spm


text:
% program demo_tutor.m
%
% P.M.T. Broersen, July 2004
%                  June 2005
% P.M.T. Broersen, November 2007

% see also demo_armasa for extensive example
%          demo_simple for basic example PSD and autocorrelation computation
%
%   Background information can be found in

%   Book:
%   Piet M.T. Broersen
%   Automatic Autocorrelation and Spectral Analysis
%   Springer-Verlag,London, 2006.
%   ISBN 1-84628-328.

%   Journal paper:
%   P. M. T. Broersen, Automatic Spectral Analysis with
%   Time Series Models, IEEE Transactions on Instrumentation
%   and Measurement, Vol. 51, No. 2, April 2002, pp. 211-216.
%
%   and many papers given in ARMASA info.txt

clc,
close all
clear all
echo on

% Generate stationary AR en MA polynomials a and b from reflection coefficients smaller than 1
% By taking reflection coefficients between -1 and +1, stationary invertible
% models are guaranteed with all poles and zeros inside the unit circle
% Reflection coefficients are transformed into parameters with rc2arset:
%

a = rc2arset([1 .3 .3])

b = rc2arset([1 -.9])

% From reflection coefficients to parameters and vice versa
[dum rc] = ar2arset(a)

N = 500;

% Generating data by starting with initial zeros is a poor idea.
% Simuarma generates N stationary data, with the asymptotical properties
% applicable from the first observation

data = simuarma(a,b,N);

% Compute and select time series models,
% with asel and bsel as parameters of the selected model
% Compute the accuracy of the selected model type and the best alternative
% model types and put them in sellog

[asel bsel sellog] = armasel(data)

% Compute the accuracy of the estimated model in comparison with
% the true process with the parameters a and b

ME = moderr(asel,bsel,a,b,N)

% Compute and plot the autocorrelation functions and spectra of the
% generating true model and the ARMAsel selected model.
% gain is the power gain or the ratio between output and input variance of ARMA process
% The length of the correlation function is chosen as 50.
%
[cortrue gain] = arma2cor(a,b,50);
corsel = arma2cor(asel,bsel,50);
[psdtrue fas] = arma2psd(a,b);
psdsel = arma2psd(asel,bsel);

echo off

figure(1),
plot(data),
xlabel('\rightarrow time axis [s]')
title([int2str(N),' ARMA(2,1) observations'])

figure(2),
plot(0:50,cortrue,0:50,corsel,'r')
title('True and estimated autocorrelation function')
xlabel('\rightarrow time lag [s]')
legend('true','estimated')

figure(3),
subplot(121),
plot(fas,psdtrue,fas,psdsel,'r')
title('True and estimated spectrum')
xlabel('\rightarrow normalized frequency \it{f/f_s}'),
legend('true','estimated',2)
ylabel('Lineair scale for PSD'),
axis tight
subplot(122),
semilogy(fas,psdtrue,fas,psdsel,'r')
title('True and estimated spectrum')
xlabel('\rightarrow normalized frequency \it{f/f_s}'),
legend('true','estimated',4)
ylabel('Log scale for PSD'),
axis tight,
subplot,

figure(4),
echo on

% The spectra on a linear scale are hardly to discern for high frequencies.
% That improves with a logaritmic scale.
% The frequency axis is between 0 and fs/2, where fs is the sampling frequency.
% Very often, fs is taken as 1 without mentioning it.
%
% Investigate the accuracy as estimated by ARMAsel for all models
% That accuracy is stored in sellog in pe_est files
% The accuracy of all models can be interpreted as the language of random data
%
% This example demonstrates how axis can be generated automatically for a
% properly scaled plot of the prediction errors
%
plot(sellog.ar.cand_order,sellog.ar.pe_est)
title('Estimated model accuracies as a function of model order and type')
xlabel('\rightarrow model order \itr')
ylabel('\rightarrow normalized model accuracy')
minar=min(sellog.ar.pe_est);
as2=axis;
as2(3)=.9*minar;
as2(4)=sellog.ar.pe_est(end);
axis(as2);
hold on,
plot(sellog.ma.cand_order,sellog.ma.pe_est,'r')
plot(sellog.arma.cand_ar_order,sellog.arma.pe_est,'g')
legend('AR(\itr\rm)','MA(\itr\rm)','ARMA(\itr,r\rm-1)',4),
hold off,

% A similar result in one line can be obtained with plotpe(sellog),
% where sellog is the third output argument of armasel

plotpe(sellog)

% Retrieve the reflection coefficients of the longest AR model that has been used in ARMAsel
% The longest AR model is allways saved as ASAglob_rc!!
%
% Recent versions of ARMASA (after 2008) save that model additionally as sellog.ar.rcarlong
%
rcarlong = ASAglobretr('ASAglob_rc');

% Determine the highest AR order that has been estimated with ARMAsel.
% Transform reflection coefficients rcarlong into a parameter vector.
% Plot spectrum (PSD) and autocorrelation function of selected and of long model.
%
% The correlation function and the PSD of the long model arlong
% are very irregular.
% Selection of the model order and type gives results with only
% statisatically significant details included.
%
maximum_order = length(rcarlong)-1
arlong = rc2arset(rcarlong);
corlang = arma2cor(arlong,1,50);
psdlang = arma2psd(arlong,1); echo off

% In arma2cor and arma2psd, 1 is used as the second parameter vector.
% That is the MA vector for true AR processes.
% The programs expect both an AR and a MA parameter vector.
% The first element of the parameter vector must always be 1.
% That can be the only parameter of the vector.

figure,
plot(0:50,cortrue,0:50,corsel,0:50,corlang)
title('True and estimated autocorrelation function')
xlabel('\rightarrow time lag [s]'),
legend('true','estimated','ARlong')

% Linear and logarithmic spectra of the same data demonstrate a strong
% preference for logarithmic plots of the PSD (power spectral density)

figure,
subplot(121),
plot(fas,psdtrue,fas,psdsel,fas,psdlang)
title('True and estimated spectrum')
xlabel('\rightarrow normalized frequency \it{f/f_s}'),
legend('true','estimated','ARlong',2)
ylabel('Linear scale for PSD'),
axis tight
subplot(122),
semilogy(fas,psdtrue,fas,psdsel,fas,psdlang)
title('True and estimated spectrum')
xlabel('\rightarrow normalized frequency \it{f/f_s}'),
legend('true','estimated','ARlong',4)
ylabel('Log scale for PSD'),
axis tight,
subplot

echo on

% Demo of variation of input variables for candidate orders
% Computation of models with prescribed type and order
%
echo off
disp('AR(10) with ar10 = armasel(data,10,0,0)')
ar10 = armasel(data,10,0,0),           % AR(10)
disp(' ')
disp('AR selected from orders between 2 and 20 with arselect= armasel(data,2:20,0,0)')
arselect = armasel(data,2:20,0,0),    % AR(selected between candidate orders 2 and 20)
disp(' ')
disp('AR(10) with ar10fromarlong = rc2arset(rcarlong,10)')
ar10fromarlong = rc2arset(rcarlong,10) % AR(10) computed with first 10 reflection coefficients
disp(' ')
disp('MA(6) model with ma6 = armasel(data,0,6,0)')
[dummy ma6] = armasel(data,0,6,0),     % MA(6)
disp(' ')
disp('ARMA(2,1) model with  armasel(data,0,0,2)')
[ar2 ma1] = armasel(data,0,0,2),       % ARMA(2,1)
disp(' ')
disp('ARMA(4,3) model with  armasel(data,0,0,4)')
[ar4 ma3] = armasel(data,0,0,4),       % ARMA(4,3)
disp(' ')
disp('Selection between candisdates AR(2), MA(3) and ARMA(5,4) with  armasel(data,2,3,5)')
[arx brx] = armasel(data,2,3,5)        % selects between the candidates AR(2), MA(3) and ARMA(5,4)
disp(' ')

% Accuracy of these models and of the selected model

disp(['ARMAsel selected from all candidates for N = ',int2str(N),' observations: ARMA(', ...
    int2str(length(asel)-1),',',int2str(length(bsel)-1),')'])
ME_selected = moderr(asel,bsel,a,b,N)
ME_ar10 = moderr(ar10,1,a,b,N)
ME_ma6 = moderr(1,ma6,a,b,N)
ME_arma21 = moderr(ar2,ma1,a,b,N)
ME_arma43 = moderr(ar4,ma3,a,b,N)
ME_arlong = moderr(arlong,1,a,b,N)

echo on,

% Prediction of the observations after the data with the selected model with arma2pred
% Lpred is the prediction horizon
% Normalizing the accuracy with the power gain
% Take care because the mean has been subtracted automatically,
% it should be added again for a correct prediction of the future
% pred_acc is the variance of the predictions
% asel and bsel have been found before with armasel in line 35
%
Lpred = 15
[pred pred_acc] = arma2pred(asel,bsel,data,Lpred);
[corsel gainsel]= arma2cor(asel,bsel)
pred_acc = pred_acc*std(data)^2/gainsel; echo off

figure,
plot(0:Lpred,[data(end) pred+mean(data)],'r*',1:Lpred, ...
    pred+1.96*sqrt(pred_acc)+mean(data),1:Lpred,pred-1.96*sqrt(pred_acc)+mean(data))
title('Prediction for \it{t > N} \rmof continuation of data \it{x_n}\rm , with 95 % confidence interval')
xlabel('\rightarrow Prediction starting with the last observation \it{x_N} \rmof the data')

echo on,

%
% Using reduced statistics
%
% REDUCED STATISTICS REDUCED STATISTICS REDUCED STATISTICS REDUCED STATISTICS
%_____________________________________________________________________________
%
% Computes the best model for the data if only a long AR model is available
% and the sample size N of the observations that were used to compute the long AR model
% Generally, there is a only small difference between ARMA models estimated from data and
% ARMA models estimated with reduced statistics.
%
%
% S. de Waele
% Automatic Spectral Analysis,
% This and more extensions of ARMASA can be found under
% spectral analysis at the download address
% http: www.mathworks.com/matlabcentral/fileexchange/

% Also programs for equidistant data where some data are missing can be
% found and programs for irregular data, all at various places in
% http: www.mathworks.com/matlabcentral/fileexchange/
% by the authors Piet Broersen and Stijn de Waele
%
% An error message can be found if the additional software is not down loaded
%
echo off

try
    disp('The model selected from armasel applied to the data was')
    asel
    bsel
    disp(' ')
    disp('    [asel_rs bsel_rs sellog_rs] = armasel_rs(arlong,N)  gives:')
    disp(' ')
    [asel_rs bsel_rs sellog_rs] = armasel_rs(arlong,N)
    disp(' ')
    disp('   The accuracy of the selected reduced statistics model')
    ME_rs = moderr(asel_rs,bsel_rs,a,b,N)
    disp(' ')
    disp('  The difference between the selected data model and the reduced statistics model')
    ME_dif = moderr(asel_rs,bsel_rs,asel,bsel,N)

    disp(' ')
    disp(' Computation of the ARMA(3,2) model with the reduced statistics algorithm')
    disp(' ')
    [a_rs b_rs sellog32_rs] = armasel_rs(arlong,N,0,0:1,3) %ARMA(3,2)
    disp(' ')
    disp('   The accuracy of the ARMA(3,2) model')
    ME_arma32_rs = moderr(a_rs,b_rs,a,b,N)

    % some care is required for giving a fixed MA or ARMA order in armasel_rs
    % using only candidate order 0 might give error messages
    % the example shows how those messages are suppressed
catch
    disp(' The program armasel_rs requires additional software of Stijn de Waele')
    disp(' Automatic Spectral Analysis')
    disp(' This can be found under spectral analysis at the download address')
    disp(' http: www.mathworks.com/matlabcentral/fileexchange/')
end


text:
function demo_excludeData ()
% // VBA toolbox //////////////////////////////////////////////////////////
%
% [posterior, out] = demo_excludeData ()
% Demo of simulation and inversion of a dynamical system that is sampled on
% an irregular grid (not all timepoints are observed)
%
% More generally, this demo shows how to deal with data exclusion
%
% /////////////////////////////////////////////////////////////////////////

% number of simulated points
N = 257;

%% Specify the model
% =========================================================================
% For the sake of the demonstration, we use a simple linear dynamical
% system $ dx/dt = Ax + phi' u $ that is directly observed $ y = x + noise $

% evolution function
% -------------------------------------------------------------------------
options.inF.A = [- 4, - 16; 4, - 4];
options.inF.dt = 1e-1;

function fx = f_evolution(x,P,u,in)
    xdot = in.A * x + diag(P) * u;
    fx = x + in.dt * xdot;
end

f_fname = @f_evolution;

% observation function
% -------------------------------------------------------------------------
g_fname = @g_Id;

%% Simulate data
% =========================================================================

% inputs
u = randn(2, N);

% parameters
x0 = [0; 0]; % initial state
theta = [1; 2]; % effect of inputs
phi = []; % no observation parameters
alpha = Inf; % deterministic
sigma = 1; % observation noise

% Build full time series of hidden states and observations
[y,x,x0,eta,e] = VBA_simulate (N,f_fname,g_fname,theta,phi,u,alpha,sigma,options,x0);

% display full time series of hidden states and observations
displaySimulations(y,x,eta,e);

%% Decimation
% =========================================================================
% here we will simulate a sampling of data on an irregular grid  by
% 'removing' data from the full simulation

% sampling grid (exponential timestps)
timepoints = 2 .^ (1 : floor (log (N) ./ log (2)));

% exclude all data but timepoints:
% if options.isYout = 1, the datapoint is ignored
options.isYout = ones (size (y));
options.isYout(:, timepoints) = 0;

%% Inversion
% =========================================================================

% dimensions of the problem
dim.n_theta = 2;
dim.n_phi = 0;
dim.n = 2;
dim.p = 2;

% Invert deterministic model
[posterior,out] = VBA_NLStateSpaceModel(y,u,f_fname,g_fname,dim,options);

% display
displayResults(posterior,out,y-e,x,x0,theta,phi,alpha,sigma);

end

author: MBB-team

repo: VBA-toolbox


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MRiLab auto generated file: DO NOT EDIT!     %
% Generated by MRiLab "DoWriteXML2m" Generator %
% MRiLab Version 1.3                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Obj, VObj]=VObj_SphereHead(flag)
%====================================
VObj.Gyro=[267538030.3797];
AttributeOpt={'Normal','MT','ME','GM'};
VObj.Model=AttributeOpt{1};
VObj.Name='VObj_SphereHead';
VObj.Notes='3D Sphere Head Object';
VObj.Type='Water';
VObj.TypeNum=[1];
VObj.XDim=[90];
VObj.XDimRes=[0.002];
VObj.YDim=[90];
VObj.YDimRes=[0.002];
VObj.ZDim=[90];
VObj.ZDimRes=[0.002];
[x,y,z]=meshgrid(((-(VObj.XDim-1)/2):((VObj.XDim-1)/2))*VObj.XDimRes,((-(VObj.YDim-1)/2):((VObj.YDim-1)/2))*VObj.YDimRes,((-(VObj.ZDim-1)/2):((VObj.ZDim-1)/2))*VObj.ZDimRes);
rho=zeros([size(x) VObj.TypeNum]);
t1=zeros([size(x) VObj.TypeNum]);
t2=zeros([size(x) VObj.TypeNum]);
t2star=zeros([size(x) VObj.TypeNum]);
chemshift=zeros(1,VObj.TypeNum);
typeflag=zeros(1,VObj.TypeNum);
lineshapeflag=zeros(1,VObj.TypeNum);
k=zeros([size(x) VObj.TypeNum^2]);
econ=zeros([size(x) 3]);
massden=zeros(size(x));
%====================================
Alpha=[1];
Color='b';
Notes='Cube Object';
p.CenterX=[0.015];
p.CenterY=[0.015];
p.CenterZ=[0.015];
p.Length=[0.03];
ChemShift=[0];
ECon=[0.5 0.5 0.5];
K=[10];
LineShapeFlag=[0];
MassDen=[1000];
Rho=[1];
T1=[1];
T2=[0.1];
T2Star=[0.01];
TypeFlag=[0];
TypeIdx=[1];
[mask, fvc]=VObjCube(p,x,y,z,flag);
Obj{1,1}=fvc;
Obj{2,1}=Color;
Obj{3,1}=Alpha;
rho(:,:,:,TypeIdx)=rho(:,:,:,TypeIdx)+mask.*Rho;
t1(:,:,:,TypeIdx)=t1(:,:,:,TypeIdx)+mask.*T1;
t2(:,:,:,TypeIdx)=t2(:,:,:,TypeIdx)+mask.*T2;
t2star(:,:,:,TypeIdx)=t2star(:,:,:,TypeIdx)+mask.*T2Star;
massden=massden+mask.*MassDen;
econ(:,:,:,1)=econ(:,:,:,1)+mask.*ECon(1);
econ(:,:,:,2)=econ(:,:,:,2)+mask.*ECon(2);
econ(:,:,:,3)=econ(:,:,:,3)+mask.*ECon(3);
chemshift(TypeIdx)=ChemShift;
typeflag(TypeIdx)=TypeFlag;
lineshapeflag(TypeIdx)=LineShapeFlag;
if ~strcmp(VObj.Model,'Normal')
for i=1:VObj.TypeNum
k(:,:,:,(TypeIdx-1)*VObj.TypeNum+i)=k(:,:,:,(TypeIdx-1)*VObj.TypeNum+i)+mask.*K(i);
end
end
p=[];
%--------------------
Alpha=[1];
Color='r';
Notes='Sphere Object';
p.CenterX=[0];
p.CenterY=[-0.04];
p.CenterZ=[0];
p.FaceNum=[20];
p.Radius=[0.04];
ChemShift=[0];
ECon=[0.5 0.5 0.5];
K=[10];
LineShapeFlag=[0];
MassDen=[1000];
Rho=[1];
T1=[1];
T2=[0.1];
T2Star=[0.01];
TypeFlag=[0];
TypeIdx=[1];
[mask, fvc]=VObjSphere(p,x,y,z,flag);
Obj{1,2}=fvc;
Obj{2,2}=Color;
Obj{3,2}=Alpha;
rho(:,:,:,TypeIdx)=rho(:,:,:,TypeIdx)+mask.*Rho;
t1(:,:,:,TypeIdx)=t1(:,:,:,TypeIdx)+mask.*T1;
t2(:,:,:,TypeIdx)=t2(:,:,:,TypeIdx)+mask.*T2;
t2star(:,:,:,TypeIdx)=t2star(:,:,:,TypeIdx)+mask.*T2Star;
massden=massden+mask.*MassDen;
econ(:,:,:,1)=econ(:,:,:,1)+mask.*ECon(1);
econ(:,:,:,2)=econ(:,:,:,2)+mask.*ECon(2);
econ(:,:,:,3)=econ(:,:,:,3)+mask.*ECon(3);
chemshift(TypeIdx)=ChemShift;
typeflag(TypeIdx)=TypeFlag;
lineshapeflag(TypeIdx)=LineShapeFlag;
if ~strcmp(VObj.Model,'Normal')
for i=1:VObj.TypeNum
k(:,:,:,(TypeIdx-1)*VObj.TypeNum+i)=k(:,:,:,(TypeIdx-1)*VObj.TypeNum+i)+mask.*K(i);
end
end
p=[];
%--------------------
VObj.Rho=rho;
VObj.T1=t1;
VObj.T2=t2;
VObj.T2Star=t2star;
VObj.ECon=econ;
VObj.MassDen=massden;
VObj.ChemShift=chemshift;
if strcmp(VObj.Model,'GM')
k=reshape(k,[size(x), VObj.TypeNum, VObj.TypeNum]);
VObj.TypeFlag=typeflag;
VObj.LineShapeFlag=lineshapeflag;
end
if ~strcmp(VObj.Model,'Normal')
VObj.K=k;
end
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function most_ex1_ed(quiet)
%MOST_EX1_ED  Examples of deterministic economic dispatch.

%   MOST
%   Copyright (c) 2015-2016, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MOST.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/most for more info.

%% set up options
define_constants;
verbose = 0;
casefile = 'ex_case3a';
mpopt = mpoption('verbose', verbose);
mpopt = mpoption(mpopt, 'out.gen', 1);
mpopt = mpoption(mpopt, 'model', 'DC');
mpopt = mpoption(mpopt, 'opf.dc.solver', 'MIPS');
mpopt = mpoption(mpopt, 'most.solver', mpopt.opf.dc.solver);
if ~verbose
    mpopt = mpoption(mpopt, 'out.all', 0);
end

%%-----  economic dispatch (no network)  -----
%% runopf
mpc = loadcase(casefile);
mpc.branch(:, RATE_A) = 0;  % disable line flow limits (mimic no network case)
r1 = rundcopf(mpc, mpopt);
Pg1 = r1.gen(:, PG);        % active generation
lam1 = r1.bus(:, LAM_P);    % nodal energy price

%% most
mpc = loadcase(casefile);
mpopt = mpoption(mpopt, 'most.dc_model', 0);    % use model with no network
mdi = loadmd(mpc);
mdo = most(mdi, mpopt);
ms = most_summary(mdo);
r2 = mdo.flow.mpc;
Pg2 = r2.gen(:, PG);        % active generation
lam2 = r2.bus(:, LAM_P);    % nodal energy price

%% comparison
Pg = [Pg1 Pg2]
lam = [lam1 lam2]

%%-----  economic dispatch (w/reserves)  -----
%% runopf_w_res
mpc = loadcase(casefile);
mpc.branch(:, RATE_A) = 0;  % disable line flow limits (mimic no network case)
r1 = runopf_w_res(mpc, mpopt);
Pg1 = r1.gen(:, PG);        % active generation
lam1 = r1.bus(:, LAM_P);    % nodal energy price
R1 = r1.reserves.R;         % reserve quantity
prc1 = r1.reserves.prc;     % reserve price

%% most
mpc = loadcase(casefile);
mpopt = mpoption(mpopt, 'most.dc_model', 0);    % use model with no network
mdi = loadmd(mpc);
mdi.FixedReserves = mpc.reserves;   % include fixed zonal reserves
mdo = most(mdi, mpopt);
ms = most_summary(mdo);
r2 = mdo.flow.mpc;
Pg2 = r2.gen(:, PG);        % active generation
lam2 = r2.bus(:, LAM_P);    % nodal energy price
R2 = r2.reserves.R;         % reserve quantity
prc2 = r2.reserves.prc;     % reserve price

%% comparison
Pg = [Pg1 Pg2]
lam = [lam1 lam2]
R = [R1 R2]
prc = [prc1 prc2]

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% ConvertBasler2pos - Convert Basler data to pos file format
%
%  USAGE
%
%    ConvertBasler2pos(filebasename,<options>)
%
%    filebasename   basename of the video file (should be filebasenmae.avi)
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'syncDatFile' name of binary file where sync signal is stored
%                   (default = filebasename_digitalin.dat)
%     'syncSampFq'  sampling freqeuncy of the sync signal (default= 20kHz)
%     'syncChan'    sync channel to read (default = 1)
%     'syncNbCh'    number of channels in the sync file (default = 1)
%     'posSampFq'   sampling frequency of the final pos file (after
%                   interpolation)
%    =========================================================================
%
% DEPENDENCIES:
%
%   bz_LoadBinary
%   Process_DetectLED

% Copyright (C) 2015 Adrien Peyrache
% edited by David Tingley, 2017
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.

warning('This function is now deprecated and needs to be replaced.')

function Process_ConvertBasler2Pos(fbasename,varargin)

syncDatFile = ['analogin.dat'];
syncSampFq  = 20000;
syncChan    = 1;
syncNbCh    = 1;
posSampFq   = 30; %in Hz

% Parse options
for i = 1:2:length(varargin),
  if ~isa(varargin{i},'char'),
    error(['Parameter ' num2str(i+3) ' is not a property (type ''help Process_ConvertBasler2Pos'' for details).']);
  end
  switch(lower(varargin{i})),
    case 'syncdatfile',
      syncDatFile = varargin{i+1};
      if ~isa(duration,'char')
        error('Incorrect value for property ''syncDatFile'' (type ''help Process_ConvertBasler2Pos'' for details).');
      end
    case 'syncsampfq',
      syncSampFq = varargin{i+1};
      if ~isa(frequency,'numeric')
        error('Incorrect value for property ''syncsampfq'' (type ''help Process_ConvertBasler2Pos'' for details).');
      end
    case 'syncchan',
      syncChan = varargin{i+1};
      if ~isa(start,'numeric')
        error('Incorrect value for property ''syncChan'' (type ''help Process_ConvertBasler2Pos'' for details).');
      end
		if start < 0, start = 0; end
    case 'syncnbch',
      syncNbCh = varargin{i+1};
      if ~isa(syncNbCh,'numeric')
        error('Incorrect value for property ''syncNbCh'' (type ''help Process_ConvertBasler2Pos'' for details).');
      end
    case 'possampfq',
      posSampFq = varargin{i+1};
      if ~isa(posSampFq,'numeric')
        error('Incorrect value for property ''posSampFq'' (type ''help Process_ConvertBasler2Pos'' for details).');
      end

    otherwise,
      error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help Process_ConvertBasler2Pos'' for details).']);
  end
end

    if ~exist([fbasename '.led'])
%         for thresh = [.05 .15 .3 .7 .9 1.5 5 10]
            thresh = .3;
        Process_DetectLED(fbasename,thresh);
%         end
    end
    pos = load([fbasename '.led']);

    fid = fopen(syncDatFile);
    dat = fread(fid,[syncNbCh inf],'uint16=>uint16');
    dat = dat(syncChan,:);

    % the line below does not work correctly with Intan analogin.dat files
%     dat = bz_LoadBinary(syncDatFile,'nchannels',syncNbCh,'channels',syncChan);

    t = (0:length(dat)-1)'/syncSampFq;

    %It's a pull-off, trigger is toward 0
%     dat = double(fastrms(fastrms(fastrms(dat,10),50),100)<500);
	dat = uint16(smooth(single(dat)))<1000;

% if bad camera pulses, use below
%     dat = smooth(dat,100);
%     dat = dat<0;

    dPos = find(diff(dat)==1);
    dNeg = find(diff(dat)==-1);

    if length(dPos) == length(dNeg)+1
        dPos = dPos(1:end-1);
    elseif length(dPos) > length(dNeg)+1 || length(dPos) < length(dNeg)
         keyboard
    end

    %Frame timing is the middle of shuter opening
    frameT = (t(dPos)+t(dNeg))/2;

    %The camera keeps on recording a few frames after software stopped
    %recording. So we skipthe last frames of the TTL
    if length(frameT)<size(pos,1);
        warning('Too many video frames!')
        keyboad
    elseif length(frameT)>size(pos,1);
        frameT = frameT(1:size(pos,1));
    end

    %We now interpolate the data at 60 Hz (or sampling fqcy specified in
    %arguments)
    recDuration = length(dat)/syncSampFq;

    timestamps = (0:1/posSampFq:recDuration-1/posSampFq)';
    pos(pos==-1) = NaN;
    newPos = interp1(frameT,pos,timestamps);

    behavior.positions.x = nanmean(pos(:,[2 4])');
    behavior.positions.y = nanmean(pos(:,[3 5])');
    behavior.position.z = [];
    dx = pos(:,3) - pos(:,5);
    dy = pos(:,2) - pos(:,4);
% 	ang = atan2(dy,dx)-angOffset;
% 	ang = mod(ang,2*pi);
    [orientation rho] = cart2pol(pos(:,3)-pos(:,5),pos(:,2)-pos(:,4));
    behavior.orientations.yaw = orientation;
    behavior.orientation.pitch = [];
    behavior.orientation.roll = [];
    behavior.rotationType = 'euler';
    behavior.description = '';
    behavior.events = [];

    warning('come up with a better head dir calculation...')
error('this function is unfinished and may not comform to the new data formatting standards...')
    behavior.timestamps = pos(:,1);
%     dlmwrite([fbasename '.pos'],[timestamps newPos],'delimiter','\t', 'precision', 32);
    save([fbasename '.behavior.mat'],'behavior');
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode