text:
%% CUBEAFEMQUADCURL quad curl equations on the unit cube
%
%   CUBEAFEMQUADCURL computes ND0-(CR-P0)-(linear ND) nonconforming approximations
%   of the quad curl equations in the unit cube. The mesh is refined
%   adaptively guided by a residual-based estimator.
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

close all; clear;

%% Set up
maxN = 5e5;
maxIt = 20;
theta = 0.3;
N = zeros(maxIt,1);
h = zeros(maxIt,1);
erru = zeros(maxIt,1);
erruL2 = zeros(maxIt,1);
errw = zeros(maxIt,1);
errphi = zeros(maxIt,1);
etaTotal = zeros(maxIt,1);
etaw = zeros(maxIt,1);
etaphi = zeros(maxIt,1);
etau = zeros(maxIt,1);

%% Generate initial mesh
H = pi;
[node,elem,HB] = cubemesh([0,H,0,H,0,H],H/2);
bdFlag = setboundary3(node,elem,'Dirichlet');
[node,elem,bdFlag,HB] = uniformrefine3(node,elem,bdFlag,HB);

%% PDE and options
pdeOrig = quadCurlDataSmooth1;
option.printlevel = 0;

%% AFEM cycles using only 1st Maxwell+ Stokes, 2nd Maxwell solved only for ref
for k = 1:maxIt
    fprintf('\n\n\nRefinement level %d\n', k);

    %% SOLVE the 1st Maxwell problem
    pdeCurl1.J = pdeOrig.quadcurlu;
    pdeCurl1.g_D = @(p) zeros(size(p,1),3);
    option.solver = 'mg';
    soln.w = Maxwellsaddle(node,elem,bdFlag,pdeCurl1,option,HB);

    %% SOLVE the Stokes problem

    % RHS of Stokes
    [curlwh, volume, curlbasis] = curlu3(node,elem,soln.w);
    pdeStokes.f = curlwh;
    pdeStokes.g_D = pdeOrig.curlu;
    option.solver = 'diag';
    option.printlevel = 1;
    solnStokes = Stokes3CRP0(node,elem,bdFlag,pdeStokes,option,HB);

    %% SOLVE the 2nd Maxwell problem (not used in the indicator for mesh refinement)

    % RHS os the 2nd Maxwell
    [elem2face,face] = dof3face(elem);
    NF = size(face,1); NT = size(elem,1);
    soln.phi = reshape(solnStokes.u,NF,3);

    phihCenter = zeros(NT,3); % at each elem center
    for j = 1:3 % each component
        phihj = soln.phi(:,j);
        phihj2elem = phihj(elem2face);
        phihCenter(:,j) = sum(phihj2elem,2)/4;
        % Crouzeix-Raviart shape function 1-3\lambda = 1/4 at center
    end

    [elem2edge,edge] = dof3edge(sort(elem,2)); % sort for edge element
    NE = size(edge,1);
    bt = zeros(NT,6);
    for j = 1:6 % 6 basis of the edge element
        bt(:,j) = dot(curlbasis(:,:,j),phihCenter,2).*volume;
    end
    fCurl = zeros(2*NE,1);
    fCurl(1:NE,:) = accumarray(elem2edge(:),bt(:),[NE 1]);

    % solve the last Maxwell problem
%     option.solver = 'direct';
    option.solver = 'mg';
    pdeCurl2.J = fCurl;
    pdeCurl2.g_D = pdeOrig.exactu;
    pdeCurl2.g = pdeOrig.g;
    soln.u = Maxwell1saddle(node,elem,bdFlag,pdeCurl2,option,HB);
    %% ESTIMATE
    option.includeU = false; % only use the indicator by w_h and phi_h
    soln.curlw = curlwh;
    N(k) = 2*length(soln.u)+3*length(soln.phi);
    [etaK, est] = estimatequadcurl(node,elem,soln,pdeOrig,option);
    errphi(k) = getL2error3(node,elem,pdeOrig.curlu,soln.phi);
    etaTotal(k) = sqrt(sum(etaK.^2));
    curluh = curlu3(node,elem,soln.u);
    [erru(k), errK] = getL2error3(node,elem,pdeOrig.curlu,curluh);
    erruL2(k) = getL2error3ND1(node,elem,pdeOrig.exactu,soln.u);
    etaw(k) = sqrt(sum((est.eta1).^2));
    etaphi(k) = sqrt(sum((est.eta2).^2));
    etau(k) = sqrt(sum((est.eta3).^2));

    %% Visualize
    figure(1);
    set(gcf, 'Position', [100 600 1000 300])
    subplot(1,3,1)
    showboundary3(node,elem,'z<pi/2','facealpha',0.5);

    subplot(1,3,2)
    errorE2V = accumarray(elem(:),repmat(errK,[4,1]),[size(node,1), 1]);
    showsolution3(node,elem,errorE2V,'z<pi/2 ','EdgeColor','k');

    subplot(1,3,3)
    etaE2V = accumarray(elem(:),repmat(etaK,[4,1]),[size(node,1), 1]);
    showsolution3(node,elem,etaE2V,'z<pi/2 ','EdgeColor','k');

    drawnow;

    %% MARK and REFINE
    if N(k) > maxN; break; end

    markedElem = mark(elem,etaK,theta);
    if k < maxIt
        [node,elem,bdFlag,HB] = bisect3(node,elem,markedElem,bdFlag,HB);
    end

end

%%
close all;
figure(1);
r1 = showrate(N,erru,10,'k-o');
r2 = showrate(N,errphi,10,'r-o');
r3 = showrate(N,etaTotal,10,'b-*');
r4 = showrate(N,erruL2,10,'color',[0.8,0.2,0.7],'marker','*');
% r4 = showrate(N,etaw,10,'color',[0.2, 0.8,0.7],'marker','*');
% r5 = showrate(N,etaphi,10,'color',[0.2, 0.8,0.7],'marker','*');
% r6 = showrate(N,etau,10,'color',[0.2, 0.8,0.7],'marker','*');
set(gca,'xscale','log','yscale','log');
T = title('Convergence');
XL = xlabel('$\#$ DoFs');
YL = ylabel('Error Magnitudes');
L = legend('$\Vert  \nabla\times(u-u_h)\Vert$', ...
    ['$N^{' num2str(r1) '}$'],...
    '$\Vert  \nabla\times u - \phi_h\Vert$', ...
    ['$N^{' num2str(r2) '}$'],...
    '$\eta(w_h,\phi_h)$', ...
    ['$N^{' num2str(r3) '}$'],...
    '$\Vert u - u_h\Vert$', ...
    ['$N^{' num2str(r4) '}$'],'LOCATION','Best');
%     '$\eta_1(w_h)$', ['$N^{' num2str(r4) '}$'],...
%     '$\eta_2(\phi_h)$', ['$N^{' num2str(r5) '}$'],...
%     '$\eta_3(u_h)$', ['$N^{' num2str(r6) '}$'],...

set([XL,YL,L],'Interpreter','latex','FontSize', 16);
set(T,'Interpreter','latex','FontSize',20);
set(gca,'TickLabelInterpreter', 'latex');
grid on;
set(gcf,'color','w','Position', [100 200 500 800])
drawnow;

%% sanity check
if N<1e4
    figure(2);
    subplot(1,2,1);
    option.scale = 2;
    option.plot = 'quiver3';
    center = (node(elem(:,1),:) + node(elem(:,2),:) + ...
        node(elem(:,3),:) + node(elem(:,4),:))/4;
    showvector3(center,pdeOrig.curlu(center),option);
    subplot(1,2,2);
    curluh = curlu3(node,elem,soln.u);
    showvector3(center,curluh,option);

    set(gcf,'color','w','Position', [500 600 750 300])
    drawnow;
end

repo: ifem


text:
function varargout = AlgorithmSelection(varargin)
% ALGORITHMSELECTION MATLAB code for AlgorithmSelection.fig
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////
%
%      ALGORITHMSELECTION, by itself, creates a new ALGORITHMSELECTION or raises the existing
%      singleton*.
%
%      H = ALGORITHMSELECTION returns the handle to a new ALGORITHMSELECTION or the handle to
%      the existing singleton*.
%
%      ALGORITHMSELECTION('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in ALGORITHMSELECTION.M with the given input arguments.
%
%      ALGORITHMSELECTION('Property','Value',...) creates a new ALGORITHMSELECTION or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before AlgorithmSelection_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to AlgorithmSelection_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help AlgorithmSelection

% Last Modified by GUIDE v2.5 13-Jun-2011 07:36:26

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @AlgorithmSelection_OpeningFcn, ...
                   'gui_OutputFcn',  @AlgorithmSelection_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before AlgorithmSelection is made visible.
function AlgorithmSelection_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to AlgorithmSelection (see VARARGIN)

% Choose default command line output for AlgorithmSelection
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

Algs = varargin{1};

SetAlgorithmLists(handles,Algs);

handles.Algorithms = Algs;
handles.InitialAlgs = Algs;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes AlgorithmSelection wait for user response (see UIRESUME)
uiwait(handles.figure1);

function SetAlgorithmLists(handles,Algorithms)

SelectedCount = 0;
AvailableCount = 0;
for i=1:length(Algorithms)
    if strcmp(Algorithms(i).Selected,'Yes')
        SelectedCount = SelectedCount + 1;
        SelectedNames{SelectedCount} = Algorithms(i).TextName;
    else
        AvailableCount = AvailableCount + 1;
        AvailableNames{AvailableCount} = Algorithms(i).TextName;
    end
end

if (AvailableCount > 0)
    set(handles.AvailableAlgorithms,'Value',1);
    set(handles.AvailableAlgorithms,'String',AvailableNames);
else
    set(handles.AvailableAlgorithms,'Value',0);
    set(handles.AvailableAlgorithms,'String',[]);
end

if (SelectedCount > 0)
    set(handles.SelectedAlgorithms,'Value',1);
    set(handles.SelectedAlgorithms,'String',SelectedNames);
else
    set(handles.SelectedAlgorithms,'Value',0);
    set(handles.SelectedAlgorithms,'String',[]);
end

if (SelectedCount > 0)
    set(handles.MoveLeft,'enable','on');
else
    set(handles.MoveLeft,'enable','off');
end
if (AvailableCount > 0)
    set(handles.MoveRight,'enable','on');
else
    set(handles.MoveRight,'enable','off');
end

% --- Outputs from this function are returned to the command line.
function varargout = AlgorithmSelection_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.Algorithms;
delete(handles.figure1);

% --- Executes on selection change in AvailableAlgorithms.
function AvailableAlgorithms_Callback(hObject, eventdata, handles)
% hObject    handle to AvailableAlgorithms (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns AvailableAlgorithms contents as cell array
%        contents{get(hObject,'Value')} returns selected item from AvailableAlgorithms

% --- Executes during object creation, after setting all properties.
function AvailableAlgorithms_CreateFcn(hObject, eventdata, handles)
% hObject    handle to AvailableAlgorithms (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on selection change in SelectedAlgorithms.
function SelectedAlgorithms_Callback(hObject, eventdata, handles)
% hObject    handle to SelectedAlgorithms (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns SelectedAlgorithms contents as cell array
%        contents{get(hObject,'Value')} returns selected item from SelectedAlgorithms

% --- Executes during object creation, after setting all properties.
function SelectedAlgorithms_CreateFcn(hObject, eventdata, handles)
% hObject    handle to SelectedAlgorithms (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in MoveRight.
function MoveRight_Callback(hObject, eventdata, handles)
% hObject    handle to MoveRight (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%get selection from Available List
Index = get(handles.AvailableAlgorithms,'Value');
Strings = get(handles.AvailableAlgorithms,'String');
Algorithm = Strings{Index};

%find Algorithm and change status
for i=1:length(handles.Algorithms)
    if strcmp(Algorithm,handles.Algorithms(i).TextName)
        handles.Algorithms(i).Selected = 'Yes';
    end
end

SetAlgorithmLists(handles,handles.Algorithms)

% Update handles structure
guidata(hObject, handles);

% --- Executes on button press in MoveLeft.
function MoveLeft_Callback(hObject, eventdata, handles)
% hObject    handle to MoveLeft (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%get selection from Available List
Index = get(handles.SelectedAlgorithms,'Value');
Strings = get(handles.SelectedAlgorithms,'String');
Algorithm = Strings{Index};

%find Algorithm and change status
for i=1:length(handles.Algorithms)
    if strcmp(Algorithm,handles.Algorithms(i).TextName)
        handles.Algorithms(i).Selected = 'No';
    end
end

SetAlgorithmLists(handles,handles.Algorithms)

% Update handles structure
guidata(hObject, handles);

% --- Executes on button press in OK.
function OK_Callback(hObject, eventdata, handles)
% hObject    handle to OK (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
uiresume(handles.figure1);

% --- Executes on button press in Cancel.
function Cancel_Callback(hObject, eventdata, handles)
% hObject    handle to Cancel (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

handles.Algorithms = handles.InitialAlgs;
guidata(hObject, handles);
uiresume(handles.figure1);

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function p = compile_bilinearslist(p);
% Precomputed list of bilinear expressions, used in when evaluating value
% which is done in heuristics code
Bilinears = find(p.variabletype==1);
BilinearsList = zeros(length(p.c),2);
for i = Bilinears
    vars = find(p.monomtable(i,:));
    BilinearsList(i,:) = vars(:)';
end
p.Bilinears = Bilinears;
p.BilinearsList = BilinearsList;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [ success ] = write_lin_svm( location, means, w, b, pos_lbl, neg_lbl)
%WRITE_LIN_SVR Summary of this function goes here
%   Detailed explanation goes here

    fileID = fopen(location, 'w');

    if(fileID ~= -1)

        % Write the regressor type 4 - linear SVM
        fwrite(fileID, 4, 'uint');

        writeMatrixBin(fileID, means, 6);
        writeMatrixBin(fileID, w, 6);
        fwrite(fileID, b, 'float64');

        fwrite(fileID, pos_lbl, 'float64');
        fwrite(fileID, neg_lbl, 'float64');

        fclose(fileID);

        success = true;

    else
        success = false;
    end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function location2D = fitCuboid(ptr, changeID)

%{
Author:
Jianxiong Xiao: http://mit.edu/jxiao/

Citation:
Please cite the following paper if you use this code in all possible ways.

@inproceedings{CuboidDetector,
 author = "Jianxiong Xiao and Bryan C. Russell and Antonio Torralba",
 title = "Localizing 3D Cuboids in Single-view Images",
 booktitle = "Advances in Neural Information Processing Systems (NIPS)",
 year = "2012",
 month = "December",
 address = "Lake Tahoe, USA"
}

Demo:
just run this function without any argument.

Function:

This program takes input of the image coordinate of the 7 corners of a
cuboid, and fit a 3D cuboid, and output the new image coordinate of the
corners.

If changeID is 1-7, it will ignore that corner, and use the rest of the
corners to estimate the result.
If changeID is 8, it will use all 7 corners to fit the cuboid.

Installation:

You need to have ceres (http://code.google.com/p/ceres-solver/) installed
in your machine (I only tested with version 1.3.0),
and compile the mex file by running compile.m

Cooridinate system:
    y
    ^
    |
    |    2-------------7
    |    |\             \
    |    | 1-------------4
    |    | |             |
    |    5 |             |
    |     \|             |
    |      3-------------6
    |
    |
    +----------------------------->x
   /
  /
 /
L
z

%}

if ~exist('ptr','var')

    subplot(2,1,1);
    imshow(imread('SUNprimitive.jpg'));

    ptr=[
        80+10 52+10  % add some error +10 to test how it gets back
        53 3
        179 55
        62 145
        152 6
        162 146
        37 90];

    edge = [1 2; 1 3; 1 4; 2 5; 3 5; 3 6; 2 7; 4 7; 4 6];
    for e=1:size(edge,1)
        hold on
        plot(ptr(edge(e,:),2),ptr(edge(e,:),1),'-r','LineWidth',3);
    end
    title('input');

end
if ~exist('changeID','var')
    changeID = 1;
end

ptr = ptr';

%% initialization

X = [-1    -1     -1     1     -1      1      1;
    1     1     -1     1     -1     -1      1;
    1    -1      1     1     -1      1     -1];
x = ptr;
N = 7;

if changeID<8
    chooseVector = true(1,7);
    chooseVector(changeID) = false;
    X = X(:,chooseVector);
    x = x(:,chooseVector);
    N = 6;
end

xi = x(1,:);
yi = x(2,:);

% Estimate Q = P*L:
A = [[zeros(4,N); X; ones(1,N); -X(1,:).*yi; -X(2,:).*yi; -X(3,:).*yi]'; ...
    [-X; -ones(1,N); zeros(4,N); X(1,:).*xi; X(2,:).*xi; X(3,:).*xi]'];
b = [yi -xi]';
Q = [A\b; 1];
Q = reshape(Q,4,3)';

[K,R,t] = decomposeP(Q); % P = K*R*[eye(3) -t].

%xxx = K*R*[X-repmat(t,1,N)]; xxx(1:2,:)./xxx([3 3],:)

f=mean(K([1 5]));
px = K(7);
py = K(8);
%cuboidH = 1;
%cuboidW = 1;

%{
        px_init = scoreSize(2)/2;
        py_init = scoreSize(1)/2;
        cuboidH = px_init/px;
        cuboidW = py_init/py;
        f = mean([K(1,1)/cuboidH, K(2,2)/cuboidW]);
        cuboidH = K(1,1)/f;
        cuboidW = K(2,2)/f;
        px = K(1,3)/cuboidH;
        py = K(2,3)/cuboidW;
        t(1) = t(1) * cuboidH;
        t(2) = t(2) * cuboidW;
%}
% K = [f 0 px; 0 f py; 0 0 1];
% xxx = K*R*[diag([cuboidH,cuboidW,1])*X-repmat(t,1,N)]; xxx(1:2,:)./xxx([3 3],:)
% sum(sum((xxx(1:2,:)./xxx([3 3],:) - ptr(:,chooseVector)).^2))

parameter = [f,px,py,t',1,1,0,0,0];

%% least square fitting

[location2D,err] = fitCuboidMex(ptr,parameter,R, int32(changeID));

%% visualization for the demo
if exist('edge','var')
    subplot(2,1,2);
    imshow(imread('SUNprimitive.jpg'));
    for e=1:size(edge,1)
        hold on
        plot(location2D(2,edge(e,:)),location2D(1,edge(e,:)),'-g','LineWidth',3);
    end
    title('result');
end

%%
if changeID<8
    location2D = location2D(:,changeID)';
end

author: jianxiongxiao

repo: ProfXkit


text:
% function [x1,x2,x3] = vec2array(X,m,flag)
%
% (c) NP, JM; SAFIR, Luebeck, 2006
%
% This function performs a basis change, nothin else
% it changes from the vector representation of a grid
% to its matlab-matrix representation
% X's dimension should be in {2,3}
% input:
%   - X     : intput gird (see getGrid for example)
%   - m     : size of grid
%   - flag  : kind of grid (one of
%               'centered', 'cell centered' (which means the same),
%               'stg', 'staggered',
%               'nodal')
%
% output:
% - x1,x2,x3 : the matlab-matrix representation of the grid
%   if X is a 2D grid (which means length(m) == 2), x3 is empty
%
function [x1,x2,x3] = vec2array(X,m,flag)
x1 = [];
x2 = [];
x3 = [];

flag = sprintf('%s-%d',flag,length(m));

switch flag,
  case {'centered-2','cell centered-2'},
    n   = prod(m);
    x1 = reshape(X(1:n)+0,m);
    x2 = reshape(X(n+1:2*n)+0,m);
  case {'stg-2','staggered-2'},
    n1 = (m(1)+1)*m(2);
    x1 = reshape(X(1:n1)+0,m(1)+1,m(2));
    n2 = m(1)*(m(2)+1);
    x2 = reshape(X(n1+(1:n2))+0,m(1),m(2)+1);
  case {'nodal-2'},
    n   = prod(m+1);
    x1 = reshape(X(1:n)+0,m+1);
    x2 = reshape(X(n+1:2*n)+0,m+1);

  case {'centered-3','cell centered-3'},
    n   = prod(m);
    x1 = reshape(X(1:n)+0,m);
    x2 = reshape(X(n+1:2*n)+0,m);
    x3 = reshape(X(2*n+1:3*n)+0,m);
  case {'stg-3','staggered-3'},
    n1 = (m(1)+1)*m(2)*m(3);
    x1 = reshape(X(1:n1)+0,m(1)+1,m(2),m(3));
    n2 = m(1)*(m(2)+1)*m(3);
    x2 = reshape(X(n1+(1:n2))+0,m(1),m(2)+1,m(3));
    n3 = m(1)*m(2)*(m(3)+1);
    x3 = reshape(X(n1+n2+(1:n3))+0,m(1),m(2),m(3)+1);
  case {'nodal-3'},
    n   = prod(m+1);
    x1 = reshape(X(1:n)+0,m+1);
    x2 = reshape(X(n+1:2*n)+0,m+1);
    x3 = reshape(X(2*n+1:3*n)+0,m+1);

  otherwise,
    error('nyi')
end;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [A,B] = hmxLowrankSub(Mh,I,J)
%+========================================================================+
%|                                                                        |
%|         OPENHMX - LIBRARY FOR H-MATRIX COMPRESSION AND ALGEBRA         |
%|           openHmx is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : hmxLowrankSub.m                               |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
%| ( === ) |   SYNOPSIS   :  Extract low-rank submatrix from H-Matrix     |
%|  `---'  |                                                              |
%+========================================================================+

%%% H-Matrix (recursion)
if (Mh.typ == 0)
    % Initialisation
    A = zeros(length(I),0,class(Mh.row{1}));
    B = zeros(0,length(J),class(Mh.row{1}));

    % Recursion
    for i = 1:4
        % Row indices
        [Ia,Ib] = ismember(I,Mh.row{i});
        Ib      = Ib(Ia);

        % Colums indices
        [Ja,Jb] = ismember(J,Mh.col{i});
        Jb      = Jb(Ja);

        % Computation for non empty indices
        if ~isempty(Ib) && ~isempty(Jb)
            % Children evaluation
            [Ai,Bi] = hmxLowrankSub(Mh.chd{i},Ib,Jb);

            % Low-rank addition for A
            if (size(Ai,2) > 0)
                tmp       = zeros(length(I),size(Ai,2),class(Mh.row{1}));
                tmp(Ia,:) = Ai;
                A         = [A , tmp];

                % Low-rank addition for B
                tmp       = zeros(size(Bi,1),length(J),class(Mh.row{1}));
                tmp(:,Ja) = Bi;
                B         = [B ; tmp];
            end
        end
    end

    % Recompression
    [A,B] = hmxQRSVD(A,B,Mh.tol);

%%% Compressed leaf
elseif (Mh.typ == 1)
    A = Mh.dat{1}(I,:);
    B = Mh.dat{2}(:,J);

%%% Full leaf
elseif (Mh.typ == 2)
    % Low-Rank conversion
    A = full(Mh.dat(I,J));
    B = eye(length(J));

    % Recompression
    [A,B] = hmxQRSVD(A,B,Mh.tol);

%%% Others
else
    error('hmxLowrankSub.m : unavailable case')
end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function spinit
% SPINIT   Add the sparse grid routines directories to Matlab path

% Author : Andreas Klimke, Universitaet Stuttgart
% Date   : October 8, 2007
% Version: 1.6

% ------------------------------------------------------------
% Sparse Grid Interpolation Toolbox
% Copyright (c) 2006 W. Andreas Klimke, Universitaet Stuttgart
% Copyright (c) 2007-2008 W. A. Klimke. All Rights Reserved.
% See LICENSE.txt for license.
% email: klimkeas@ians.uni-stuttgart.de
% web  : http://www.ians.uni-stuttgart.de/spinterp
% ------------------------------------------------------------

try
	disp('------------------------------------------------------------------------------');
	disp('Sparse Grid Interpolation Toolbox');
	disp('------------------------------------------------------------------------------');
	disp('Version 5.1.1');
	disp('(c) 2006 W. Andreas Klimke, Universitaet Stuttgart.');
	disp('(c) 2007-2008 W. A. Klimke. All Rights Reserved.');
	disp('Please see LICENSE.txt for license information.');
	disp('------------------------------------------------------------------------------');
	fullpath = mfilename('fullpath');
	pathstr = fileparts(fullpath);
	disp('Adding sparse grid routines directory to Matlab search path.');
	addpath(pathstr);
	disp('Adding examples directory to Matlab search path.');
	demopath = [pathstr filesep 'examples'];
	addpath(demopath);
	disp('------------------------------------------------------------------------------');
	disp(['Type >>help spinterp<< or >>help spvals<< for brief usage' ...
				' information.']);
  disp(['Type >>doc spinterptool<< to open the documentation in the' ...
        ' help browser.']);
	disp(['Type >>demo toolbox ''Sparse Grid Interpolation''<< to' ...
	      ' select demo files.']);
	disp(' ');
catch
	disp('Initialization failed. The following error occurred:');
	rethrow(lasterror);
end

author: johannesgerer

repo: jburkardt-m


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

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

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

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

predictions =  sigmoid(X*theta);

leftPart = -y' * log(predictions);

rightPart = (1 - y') * log(1 - predictions);

J = (1 / m) * (leftPart - rightPart);

grad = (1 / m) * ((predictions - y)' * X);

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

end

author: vugsus

repo: coursera-machine-learning


text:
function pf=lpcar2pf(ar,np)
%LPCAR2PF Convert AR coefs to power spectrum PF=(AR,NP)
% the power spectrum will contain np+2 points ranging from Dc to nyquist frequency
% The routine is faster if NP+1 is a power of 2.

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: lpcar2pf.m,v 1.4 2007/05/04 07:01:38 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,p1]=size(ar);
if nargin<2 np=p1-1; end
pf=abs(rfft(ar.',2*np+2).').^(-2);

author: decouples


text:
function x = spm_invXcdf(F,v)
% Inverse Cumulative Distribution Function (CDF) of Chi-squared distribution
% FORMAT x = spm_invXcdf(F,v)
%
% F - CDF (lower tail p-value)
% v - degrees of freedom (v>0, non-integer d.f. accepted)
% x - Chi-squared ordinates at which CDF F(x)=F
%__________________________________________________________________________
%
% spm_invXcdf implements the inverse Cumulative Distribution of the
% Chi-squared distribution.
%
% Definition:
%--------------------------------------------------------------------------
% The Chi-squared distribution with v degrees of freedom is defined for
% positive integer v and x in [0,Inf). The Cumulative Distribution
% Function (CDF) F(x) is the probability that a realisation of a
% Chi-squared random variable X has value less than x. F(x)=Pr{X<x}:
% (See Evans et al., Ch8)
%
% Variate relationships: (Evans et al., Ch8 & Ch18)
%--------------------------------------------------------------------------
% The Chi-squared distribution with v degrees of freedom is equivalent to
% the Gamma distribution with scale parameter 1/2 and shape parameter v/2.
%
% Algorithm:
%--------------------------------------------------------------------------
% Using routine spm_invGcdf for Gamma distribution, with appropriate
% parameters.
%
% References:
%--------------------------------------------------------------------------
% Evans M, Hastings N, Peacock B (1993)
%       "Statistical Distributions"
%        2nd Ed. Wiley, New York
%
% Abramowitz M, Stegun IA, (1964)
%       "Handbook of Mathematical Functions"
%        US Government Printing Office
%
% Press WH, Teukolsky SA, Vetterling AT, Flannery BP (1992)
%       "Numerical Recipes in C"
%        Cambridge
%
%__________________________________________________________________________
% Copyright (C) 1993-2011 Wellcome Trust Centre for Neuroimaging

% Andrew Holmes
% $Id: spm_invXcdf.m 4182 2011-02-01 12:29:09Z guillaume $

%-Check enough arguments
%--------------------------------------------------------------------------
if nargin<2, error('Insufficient arguments'), end

%-Computation
%--------------------------------------------------------------------------
x = spm_invGcdf(F,v/2,1/2);

author: spm


text:
function imo = cnn_shape_get_batch(images, varargin)
% Modified from CNN_IMAGENET_GET_BATCH
%
% - added `pad` option
% - deals with images of types other than jpeg
% - augmentation made consistent across views

opts.imageSize = [227, 227] ;
opts.border = [29, 29] ;
opts.pad = 0;  % [TOP BOTTOM LEFT RIGHT]
opts.nViews = 1;
opts.keepAspect = true ;
opts.numAugments = 1 ;
opts.transformation = 'none' ;
opts.averageImage = [] ;
opts.rgbVariance = zeros(0,3,'single') ;
opts.interpolation = 'bilinear' ;
opts.numThreads = 1 ;
opts.prefetch = false ;
opts = vl_argparse(opts, varargin);

% if only one value is given, apply the same amount of padding to all borders
if numel(opts.pad)==1, opts.pad = repmat(opts.pad,[1 4]); end
if numel(opts.border)==1, opts.border = repmat(opts.border,[1 2]); end

% fetch is true if images is a list of filenames (instead of
% a cell array of images)
fetch = numel(images) >= 1 && ischar(images{1}) ;

% isjpg is true if all images to fetch are of jpeg format
isjpg = fetch && strcmpi(images{1}(end-3:end),'.jpg');

assert(mod(numel(images),opts.nViews)==0, '''nViews'' is incompatible with input');
nViews = opts.nViews;
nShapes = numel(images)/nViews;

if opts.prefetch
  if isjpg, vl_imreadjpeg(images, 'numThreads', opts.numThreads, 'prefetch'); end
  imo = [] ;
  return ;
end

if fetch
  if isjpg,
    im = vl_imreadjpeg(images,'numThreads', opts.numThreads) ;
  else
    im = cell(size(images));
  end
else
  im = images ;
end

tfs = [] ;
switch opts.transformation
  case 'none'
    tfs = [
      .5 ;
      .5 ;
       0 ] ;
  case 'f5'
    tfs = [...
      .5 0 0 1 1 .5 0 0 1 1 ;
      .5 0 1 0 1 .5 0 1 0 1 ;
       0 0 0 0 0  1 1 1 1 1] ;
  case 'f25'
    [tx,ty] = meshgrid(linspace(0,1,5)) ;
    tfs = [tx(:)' ; ty(:)' ; zeros(1,numel(tx))] ;
    tfs_ = tfs ;
    tfs_(3,:) = 1 ;
    tfs = [tfs,tfs_] ;
  case 'stretch'
  otherwise
    error('Uknown transformations %s', opts.transformation) ;
end
[~,transformations] = sort(rand(size(tfs,2), numel(images)), 1) ;

if ~isempty(opts.rgbVariance) && isempty(opts.averageImage)
  opts.averageImage = zeros(1,1,3) ;
end
if numel(opts.averageImage) == 3
  opts.averageImage = reshape(opts.averageImage, 1,1,3) ;
end

imo = zeros(opts.imageSize(1), opts.imageSize(2), 3, ...
            numel(images)*opts.numAugments, 'single') ;

for i=1:nShapes,
  for j=1:nViews,
    % acquire image
    idx = (i-1)*nViews + j;
    if isempty(im{idx})
      imt = imread(images{idx}) ;
      imt = single(imt) ; % faster than im2single (and multiplies by 255)
    else
      imt = im{idx} ;
    end
    if size(imt,3) == 1
      imt = cat(3, imt, imt, imt) ;
    end
    if j==1,
      imArr = zeros(size(imt,1),size(imt,2),3,nViews,'single');
    end
    imArr(:,:,:,j) = imt;
  end

  % pad
  if ~isempty(opts.pad) && any(opts.pad>0),
    w = size(imArr,2);
    h = size(imArr,1);
    imArrTmp = imArr;
    imArr = 255*ones(h+sum(opts.pad(1:2)), w+sum(opts.pad(3:4)), 3, nViews, 'single');
    imArr(opts.pad(1)+(1:h), opts.pad(3)+(1:w),:,:) = imArrTmp;
  end

  % resize
  w = size(imArr,2) ;
  h = size(imArr,1) ;
  factor = [(opts.imageSize(1)+opts.border(1))/h ...
            (opts.imageSize(2)+opts.border(2))/w];

  if opts.keepAspect
    factor = max(factor) ;
  end
  if any(abs(factor - 1) > 0.0001)
    imArr = imresize(imArr, ...
                     'scale', factor, ...
                     'method', opts.interpolation) ;
  end

  % crop & flip
  w = size(imArr,2) ;
  h = size(imArr,1) ;
  for ai = 1:opts.numAugments
    switch opts.transformation
      case 'stretch'
        sz = round(min(opts.imageSize(1:2)' .* (1-0.1+0.2*rand(2,1)), [h;w])) ;
        dx = randi(w - sz(2) + 1, 1) ;
        dy = randi(h - sz(1) + 1, 1) ;
        flip = rand > 0.5 ;
      otherwise
        tf = tfs(:, transformations(mod(ai-1, numel(transformations)) + 1)) ;
        sz = opts.imageSize(1:2) ;
        dx = floor((w - sz(2)) * tf(2)) + 1 ;
        dy = floor((h - sz(1)) * tf(1)) + 1 ;
        flip = tf(3) ;
    end
    sx = round(linspace(dx, sz(2)+dx-1, opts.imageSize(2))) ;
    sy = round(linspace(dy, sz(1)+dy-1, opts.imageSize(1))) ;
    if flip, sx = fliplr(sx) ; end

    if ~isempty(opts.averageImage)
      offset = opts.averageImage ;
      if ~isempty(opts.rgbVariance)
        offset = bsxfun(@plus, offset, reshape(opts.rgbVariance * randn(3,1), 1,1,3)) ;
      end
      imo(:,:,:,(ai-1)*numel(images)+(i-1)*nViews+(1:nViews)) = bsxfun(@minus, imArr(sy,sx,:,:), offset) ;
    else
      imo(:,:,:,(ai-1)*numel(images)+(i-1)*nViews+(1:nViews)) = imArr(sy,sx,:,:) ;
    end
  end
end

author: suhangpro

repo: mvcnn

save_path: github-repos/MATLAB/suhangpro-mvcnn


text:
function y = not(x)
%NOT Logical NOT (~) for sptensors.
%
%   ~X performs a logical not on the input tensor X. The result always
%   returned as a sparse tensor.
%
%   See also SPTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Observations for sparse matrix case.
% The result of ~a is sparse.

%% Then compute those indicies that are not in x
subs = setdiff(allsubs(x),x.subs,'rows');

%% Assemble final result
y = sptensor(subs,true,x.size);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [c]=kron2(a,b)
% Computes the Pseudokronecker product of TT tensors:
%   [C]=KRON2(A,B)
%   c{i}=kron(a{i},b{i}) <- "a" indices vary faster!
%   mode sizes and ranks are multiplied
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

if ( isempty(a) )
  c=b;
  return
elseif ( isempty(b) )
  c=a;
  return
end
if ( size(a.core,1) == 1 )
   a.core=(a.core).';
end
if ( size(b.core,1) == 1 )
   b.core=(b.core).';
end

d = a.d;
if (b.d~=d)
    error('Could not compute kron2 from the tensors with not equal dimensions\n');
end;
cra = a.core;
na = a.n;
ra = a.r;
psa = a.ps;
crb = b.core;
nb = b.n;
rb = b.r;
psb = b.ps;

c=tt_tensor;
n = [(na).*(nb)];
r = [(ra).*(rb)];
ps = cumsum([1; n.*r(2:d+1).*r(1:d)]);

c.d=d;
c.n=n;
c.r=r;
c.ps = ps;
cr = zeros(ps(d+1)-1, 1);
for i=1:d
    a1 = cra(psa(i):psa(i+1)-1);
    b1 = crb(psb(i):psb(i+1)-1);
    c1 = a1*(b1.'); % Size ra1*na*ra2, rb1*nb*rb2
    c1 = reshape(c1, ra(i), na(i), ra(i+1), rb(i), nb(i), rb(i+1));
    c1 = permute(c1, [1, 4, 2, 5, 3, 6]); % We need rc1, nc*rc2
    cr(ps(i):ps(i+1)-1) = c1(:);
end;

c.core=cr;

end

author: oseledets

repo: TT-Toolbox


text:
function [ output_meta ] = meta2sicd_csm( HDF5_fid )
%META2SICD_CSM Converts CSM HDF5 into a SICD-style metadata structure
%
% Also handles KOMPSAT-5, which is nearly identical format
%
% Takes as input a file identifier to an open HDF5 file-- as is returned by
% H5F.open
%
% SICD fields not currently computed:
%    ValidData
%
% Written by: Wade Schwartzkopf, NGA/Research
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

SECONDS_IN_A_DAY = 24*60*60;

mission_id = deblank(get_hdf_attribute(HDF5_fid,'Mission ID')');

%% Setup HDF ids
numbands = H5G.get_num_objs(HDF5_fid);
if strcmp( mission_id, 'CSG')
    numbands = 1;  % CSG only has one SICD per HDF5, even for quadpol
end
for i=1:numbands % "pingpong" mode has multiple polarizations
    groupname=['/S0' num2str(i)];
    group_id(i)=H5G.open(HDF5_fid,groupname);
    switch mission_id
        case {'CSK','KMPS'}
            dset_id(i)=H5D.open(HDF5_fid,[groupname '/SBI']);
        otherwise  % 'CSG'
            dset_id(i)=H5D.open(HDF5_fid,[groupname '/IMG']);
    end
    dspace_id(i)=H5D.get_space(dset_id(i));
end

%% CollectionInfo
output_meta.CollectionInfo.CollectorName=deblank(get_hdf_attribute(HDF5_fid,'Satellite ID')');
[collectStart, collectStartFrac]=datenum_w_frac(...
    get_hdf_attribute(HDF5_fid,'Scene Sensing Start UTC')');
% output_meta.CollectionInfo.CoreName=num2str(get_hdf_attribute(HDF5_fid,'Programmed Image ID'));
output_meta.CollectionInfo.CoreName=[datestr(collectStart, 'ddmmmyy'), ...
    '_', output_meta.CollectionInfo.CollectorName, '_', ...
    datestr(collectStart, 'HHMMSS')];
output_meta.CollectionInfo.CollectType='MONOSTATIC';
acq_mode = deblank(get_hdf_attribute(HDF5_fid,'Acquisition Mode')');
switch mission_id
    case 'CSK'
        switch acq_mode
            case {'HIMAGE','PINGPONG'} % Stripmap
                output_meta.CollectionInfo.RadarMode.ModeType='STRIPMAP';
            case {'WIDEREGION','HUGEREGION'} % ScanSAR
                output_meta.CollectionInfo.RadarMode.ModeType='STRIPMAP';
            case {'ENHANCED SPOTLIGHT','SMART'} % "Spotlight"
                output_meta.CollectionInfo.RadarMode.ModeType='DYNAMIC STRIPMAP';
        end
    case 'KMPS'
        switch acq_mode
            case {'STANDARD','ENHANCED STANDARD'} % Stripmap
                output_meta.CollectionInfo.RadarMode.ModeType='STRIPMAP';
            case {'WIDE SWATH','ENHANCED WIDE SWATH'} % ScanSAR
                output_meta.CollectionInfo.RadarMode.ModeType='STRIPMAP';
            case {'HIGH RESOLUTION', 'ENHANCED HIGH RESOLUTION', 'ULTRA HIGH RESOLUTION'} % "Spotlight"
                output_meta.CollectionInfo.RadarMode.ModeType='DYNAMIC STRIPMAP';
        end
    case 'CSG'
        if strncmp(acq_mode, 'SPOTLIGHT', 9)
            output_meta.CollectionInfo.RadarMode.ModeType='DYNAMIC STRIPMAP';
        elseif any(strcmp(acq_mode, {'STRIPMAP', 'QUADPOL'}))
            output_meta.CollectionInfo.RadarMode.ModeType='STRIPMAP';
        else
            output_meta.CollectionInfo.RadarMode.ModeType='DYNAMIC STRIPMAP';
        end
end
output_meta.CollectionInfo.RadarMode.ModeID=deblank(get_hdf_attribute(HDF5_fid,'Multi-Beam ID')');
output_meta.CollectionInfo.Classification='UNCLASSIFIED';

%% ImageCreation
output_meta.ImageCreation.DateTime=datenum(get_hdf_attribute(...
    HDF5_fid,'Product Generation UTC')','yyyy-mm-dd HH:MM:SS.FFF');
output_meta.ImageCreation.Profile='Prototype';

%% ImageData
output_meta.ImageData=struct(); % Just a placeholder
% Most subfields added below in "per band" section
% Used for computing SCP later

%% GeoData
if strcmp(deblank(get_hdf_attribute(HDF5_fid,'Ellipsoid Designator')'),'WGS84')
    output_meta.GeoData.EarthModel='WGS_84';
end
% Most subfields added below in "per band" section

%% Grid
if strncmp(get_hdf_attribute(HDF5_fid,'Projection ID')','SLANT',5)
    output_meta.Grid.ImagePlane='SLANT';
    output_meta.Grid.Type='RGZERO';
else
    output_meta.Grid.ImagePlane='GROUND';
end
output_meta.Grid.Row.Sgn=-1; % Always true for CSM
output_meta.Grid.Col.Sgn=-1; % Always true for CSM
fc=get_hdf_attribute(HDF5_fid,'Radar Frequency'); % Center frequency
output_meta.Grid.Row.KCtr=2*fc/SPEED_OF_LIGHT;
output_meta.Grid.Col.KCtr=0;
output_meta.Grid.Row.DeltaKCOAPoly=0;
output_meta.Grid.Row.WgtType.WindowName=upper(deblank(get_hdf_attribute(HDF5_fid,'Range Focusing Weighting Function')'));
if strcmpi(output_meta.Grid.Row.WgtType.WindowName,'HAMMING') % The usual CSM weigting
    rowval = get_hdf_attribute(HDF5_fid,'Range Focusing Weighting Coefficient');
elseif strcmpi(output_meta.Grid.Row.WgtType.WindowName,'GENERAL_COSINE')
    output_meta.Grid.Row.WgtType.WindowName='HAMMING';  % They are effectively the same
    rowval = 1-get_hdf_attribute(HDF5_fid,'Range Focusing Weighting Coefficient');
end
if exist('rowval','var')
    output_meta.Grid.Row.WgtType.Parameter.name = 'COEFFICIENT';
    output_meta.Grid.Row.WgtType.Parameter.value = num2str(rowval);
    output_meta.Grid.Row.WgtFunct = raised_cos_fun(512,rowval);
end
output_meta.Grid.Col.WgtType.WindowName=upper(deblank(get_hdf_attribute(HDF5_fid,'Azimuth Focusing Weighting Function')'));
if strcmpi(output_meta.Grid.Col.WgtType.WindowName,'HAMMING') % The usual CSM weigting
    colval = get_hdf_attribute(HDF5_fid,'Azimuth Focusing Weighting Coefficient');
elseif strcmpi(output_meta.Grid.Col.WgtType.WindowName,'GENERAL_COSINE')
    output_meta.Grid.Col.WgtType.WindowName='HAMMING';  % They are effectively the same
    colval = 1-get_hdf_attribute(HDF5_fid,'Azimuth Focusing Weighting Coefficient');
end
if exist('colval','var')
    output_meta.Grid.Col.WgtType.Parameter.name = 'COEFFICIENT';
    output_meta.Grid.Col.WgtType.Parameter.value = num2str(colval);
    output_meta.Grid.Col.WgtFunct = raised_cos_fun(512,colval);
end
% More subfields added below in "per band" section

%% Timeline
[collectStart, collectStartFrac]=datenum_w_frac(...
    get_hdf_attribute(HDF5_fid,'Scene Sensing Start UTC')');
[collectEnd, collectEndFrac]=datenum_w_frac(...
    get_hdf_attribute(HDF5_fid,'Scene Sensing Stop UTC')');
% We loose a bit of precision when assigning the SICD CollectStart
% field, since a MATLAB serial date number just doesn't have enough
% bits to handle the full precision given in the TSX XML. However, all
% relative times within the SICD metadata structure will be computed at
% full precision.
output_meta.Timeline.CollectStart=collectStart + (collectStartFrac/SECONDS_IN_A_DAY);
output_meta.Timeline.CollectDuration=...
    round((collectEnd-collectStart)*SECONDS_IN_A_DAY) + ... % Convert days to seconds
    (collectEndFrac-collectStartFrac); % Handle fractional seconds
output_meta.Timeline.IPP.Set.TStart=0;
output_meta.Timeline.IPP.Set.TEnd=0; % Apply real value later.  Just a placeholder.
output_meta.Timeline.IPP.Set.IPPStart=uint32(0);
% More subfields added below in "per band" section

%% Position
% Compute polynomial from state vectors
[ref_time, ref_time_frac]=datenum_w_frac(...
    get_hdf_attribute(HDF5_fid,'Reference UTC')');
% Times in SICD are with respect to time from start of collect, but
% time in CSM are generally with respect to reference time.
ref_time_offset = round((ref_time-collectStart)*SECONDS_IN_A_DAY) + ... % Convert from days to secs
    (ref_time_frac-collectStartFrac); % Handle fractional seconds
state_vector_T  = get_hdf_attribute(HDF5_fid,'State Vectors Times').'; % In seconds
state_vector_T = state_vector_T + ref_time_offset; % Make with respect to Timeline.CollectStart
state_vector_pos  = get_hdf_attribute(HDF5_fid,'ECEF Satellite Position');
% sv2poly.m shows ways to determine best polynomial order, but 5th is almost always best
polyorder=min(5, numel(state_vector_T) - 1);
P_x  = polyfit(state_vector_T, state_vector_pos(1,:), polyorder);
P_y  = polyfit(state_vector_T, state_vector_pos(2,:), polyorder);
P_z  = polyfit(state_vector_T, state_vector_pos(3,:), polyorder);
% We don't use these since they are derivable from the position polynomial
% state_vector_vel = get_hdf_attribute(HDF5_fid,'ECEF Satellite Velocity');
% state_vector_acc = get_hdf_attribute(HDF5_fid,'ECEF Satellite Acceleration');
% P_vx = polyfit(state_vector_T, state_vector_vel(1,:), polyorder);
% P_vy = polyfit(state_vector_T, state_vector_vel(2,:), polyorder);
% P_vz = polyfit(state_vector_T, state_vector_vel(3,:), polyorder);
% P_ax = polyfit(state_vector_T, state_vector_acc(1,:), polyorder);
% P_ay = polyfit(state_vector_T, state_vector_acc(2,:), polyorder);
% P_az = polyfit(state_vector_T, state_vector_acc(3,:), polyorder);
% Store position polynomial
output_meta.Position.ARPPoly.X  = P_x(end:-1:1).';
output_meta.Position.ARPPoly.Y  = P_y(end:-1:1).';
output_meta.Position.ARPPoly.Z  = P_z(end:-1:1).';

%% RadarCollection
switch mission_id
    case {'CSK', 'KMPS'}
        pols=cell(numbands,1);
        for i=1:numbands
            pols{i}=deblank(get_hdf_attribute(group_id(i),'Polarisation')');
        end
    case 'CSG'
        if strcmp(acq_mode, 'QUADPOL')
            pols = {'HH', 'HV', 'VH', 'VV'};
        else
            pols = {deblank(get_hdf_attribute(HDF5_fid,'Polarization')')};
        end
end
for i = 1:numel(pols)
    output_meta.RadarCollection.RcvChannels.ChanParameters(i).TxRcvPolarization=[pols{i}(1) ':' pols{i}(2)];
end
tx_pol = [pols{:}]; tx_pol = unique(tx_pol(1:2:end));
if numel(tx_pol)==1
    output_meta.RadarCollection.TxPolarization = tx_pol;
else
    output_meta.RadarCollection.TxPolarization = 'SEQUENCE';
    for i = 1:numel(tx_pol)
        output_meta.RadarCollection.TxSequence(i).TxStep = i;
        output_meta.RadarCollection.TxSequence(i).TxPolarization = tx_pol(i);
    end
end
% Most subfields added below in "per band" section

%% ImageFormation
output_meta.ImageFormation.RcvChanProc=struct('NumChanProc',1,'PRFScaleFactor',1);
output_meta.ImageFormation.ImageFormAlgo='RMA';
output_meta.ImageFormation.TStartProc=0;
output_meta.ImageFormation.TEndProc=output_meta.Timeline.CollectDuration;
output_meta.ImageFormation.STBeamComp='NO';
output_meta.ImageFormation.ImageBeamComp='SV';
output_meta.ImageFormation.AzAutofocus='NO';
output_meta.ImageFormation.RgAutofocus='NO';
% More subfields added below in "per band" section
output_meta.RMA.RMAlgoType='OMEGA_K';
output_meta.RMA.ImageType='INCA';
output_meta.RMA.INCA.FreqZero=fc;

%% SCPCOA
output_meta.SCPCOA.SideOfTrack=get_hdf_attribute(HDF5_fid,'Look Side')';
output_meta.SCPCOA.SideOfTrack=upper(output_meta.SCPCOA.SideOfTrack(1));
% Most subfields added below in "per band" section, but we grab this field
% now so we know to flip from CSM's EARLY-LATE column order to SICD's
% view-from-above column order.

%% Process fields specific to each polarimetric band
band_independent_meta=output_meta; % Values that are consistent across all bands
grouped_meta=cell(numbands,1);
for i=1:numbands
    output_meta=band_independent_meta;

    %% Doppler polynomials
    % These polynomials are used later to determine RMA.INCA.DopCentroidPoly
    try
        dop_rate_poly_rg=get_hdf_attribute(HDF5_fid,'Doppler Rate vs Range Time Polynomial');
    catch
        dop_rate_poly_rg=get_hdf_attribute(group_id(i),'Doppler Rate vs Range Time Polynomial');
    end
    dop_rate_poly_rg=dop_rate_poly_rg(1:find(dop_rate_poly_rg~=0,1,'last'));  % strip off zero coefficients
    dop_rate_poly_rg=-sign(dop_rate_poly_rg(1))*dop_rate_poly_rg;  % Must be negative (KS5 makes it positive)
    switch mission_id
        case {'CSK', 'KMPS'}
            t_az_ref = get_hdf_attribute(HDF5_fid,'Azimuth Polynomial Reference Time');
            t_rg_ref = get_hdf_attribute(HDF5_fid,'Range Polynomial Reference Time');
            dop_poly_az=get_hdf_attribute(HDF5_fid,'Centroid vs Azimuth Time Polynomial');
            dop_poly_rg=get_hdf_attribute(HDF5_fid,'Centroid vs Range Time Polynomial');
            dop_poly_rg=dop_poly_rg(1:find(dop_poly_rg~=0,1,'last'));  % strip off zero coefficients
        case 'CSG'
            B0001_meta=H5G.open(group_id(i),'B0001');
            az_ref_time_nozd = get_hdf_attribute(group_id(i),'Azimuth Polynomial Reference Time');
            first_time = get_hdf_attribute(B0001_meta,'Azimuth First Time');
            last_time = get_hdf_attribute(B0001_meta,'Azimuth Last Time');
            az_fit_times = linspace(first_time, last_time, 11);

            geom_dop_cent_poly = get_hdf_attribute(group_id(i),'Doppler Centroid vs Azimuth Time Polynomial - RAW');

            try
                dop_rate_poly = get_hdf_attribute(HDF5_fid,'Doppler Rate vs Azimuth Time Polynomial');
            catch
                dop_rate_poly = get_hdf_attribute(group_id(i),'Doppler Rate vs Azimuth Time Polynomial');
            end

            centroid_values = polyval(geom_dop_cent_poly(end:-1:1), az_fit_times - az_ref_time_nozd);
            rate_values = polyval(dop_rate_poly(end:-1:1), az_fit_times - az_ref_time_nozd);
            zd_times = az_fit_times - centroid_values ./ rate_values;
            t_az_ref = get_hdf_attribute(group_id(i),'Azimuth Polynomial Reference Time - ZD');
            dop_poly_az = polyfit(zd_times - t_az_ref, centroid_values, 4);
            dop_poly_az = dop_poly_az(end:-1:1)';
            dop_poly_az = dop_poly_az(1:find(dop_poly_az~=0,1,'last')); % strip off zero coefficients

            t_rg_ref = get_hdf_attribute(group_id(i),'Range Polynomial Reference Time');
            dop_poly_rg = get_hdf_attribute(group_id(i),'Doppler Centroid vs Range Time Polynomial');
            dop_poly_rg = dop_poly_rg(1:find(dop_poly_rg~=0,1,'last'));  % strip off zero coefficients
    end

    %% ImageData
    [num_dims datasize] = H5S.get_simple_extent_dims(dspace_id(i)); % All polarizations should be same size
    output_meta.ImageData.NumCols=uint32(datasize(1));
    output_meta.ImageData.NumRows=uint32(datasize(2));
    output_meta.ImageData.FullImage=output_meta.ImageData;
    output_meta.ImageData.FirstRow=uint32(0); output_meta.ImageData.FirstCol=uint32(0);
    switch H5T.get_class(H5D.get_type(dset_id(i)))
        case 0  % integer
            output_meta.ImageData.PixelType='RE16I_IM16I';
        case 1  % float
            output_meta.ImageData.PixelType='RE32F_IM32F';
    end
    % There are many different options for picking the SCP point.  We chose
    % the point that is closest to the reference zero-doppler and range
    % times in the CSM metadata.
    t_az_first=get_hdf_attribute(dset_id(i),'Zero Doppler Azimuth First Time'); % Zero doppler time of first column
    ss_az_s=get_hdf_attribute(dset_id(i),'Line Time Interval'); % Image column spacing in zero doppler time (seconds)
    output_meta.ImageData.SCPPixel.Col=uint32(round((t_az_ref-t_az_first)/ss_az_s) + 1);
    if output_meta.SCPCOA.SideOfTrack=='L'
        % Order of columns in SICD goes in reverse time for left-looking
        ss_az_s=-ss_az_s;
        output_meta.ImageData.SCPPixel.Col=uint32(output_meta.ImageData.NumCols-output_meta.ImageData.SCPPixel.Col-1);
        % First column in SICD is actually last line in CSM terminology
        t_az_first=get_hdf_attribute(dset_id(i),'Zero Doppler Azimuth Last Time');
    end
    t_rg_first=get_hdf_attribute(dset_id(i),'Zero Doppler Range First Time'); % Range time of first row
    if strncmp(get_hdf_attribute(HDF5_fid,'Product Type'),'SCS',3)
        % 'Column Time Interval' does not exist in detected products.
        ss_rg_s=get_hdf_attribute(dset_id(i),'Column Time Interval'); % Row spacing in range time (seconds)
        output_meta.ImageData.SCPPixel.Row=uint32(round((t_rg_ref-t_rg_first)/ss_rg_s) + 1);
    end
    % How Lockheed seems to pick the SCP:
    output_meta.ImageData.SCPPixel.Col = floor(datasize(1)/2);
    output_meta.ImageData.SCPPixel.Row = ceil(datasize(2)/2)-1;

    %% GeoData
    % Initially, we just seed this with a rough value.  Later we will put
    % in something more precise.
    switch mission_id
        case {'CSK', 'KMPS'}
            latlon=get_hdf_attribute(group_id(i),'Centre Geodetic Coordinates');
        case 'CSG'
            latlon=get_hdf_attribute(HDF5_fid,'Scene Centre Geodetic Coordinates');
    end
    output_meta.GeoData.SCP.LLH.Lat=latlon(1);
    output_meta.GeoData.SCP.LLH.Lon=latlon(2);
    % CSM generally gives HAE as zero.  Perhaps we should adjust this to DEM.
    output_meta.GeoData.SCP.LLH.HAE=latlon(3);
    ecf=geodetic_to_ecf([output_meta.GeoData.SCP.LLH.Lat output_meta.GeoData.SCP.LLH.Lon output_meta.GeoData.SCP.LLH.HAE]);
    output_meta.GeoData.SCP.ECF.X=ecf(1);
    output_meta.GeoData.SCP.ECF.Y=ecf(2);
    output_meta.GeoData.SCP.ECF.Z=ecf(3);
    % Calling derived_sicd_fields at the end will populate these fields
    % with the sensor model, so we don't need to do it here.
    % latlon=get_hdf_attribute(dset_id(i),'Top Left Geodetic Coordinates');
    % output_meta.GeoData.ImageCorners.ICP.FRFC.Lat=latlon(1);
    % output_meta.GeoData.ImageCorners.ICP.FRFC.Lon=latlon(2);
    % latlon=get_hdf_attribute(dset_id(i),'Bottom Left Geodetic Coordinates');
    % output_meta.GeoData.ImageCorners.ICP.FRLC.Lat=latlon(1);
    % output_meta.GeoData.ImageCorners.ICP.FRLC.Lon=latlon(2);
    % latlon=get_hdf_attribute(dset_id(i),'Bottom Right Geodetic Coordinates');
    % output_meta.GeoData.ImageCorners.ICP.LRLC.Lat=latlon(1);
    % output_meta.GeoData.ImageCorners.ICP.LRLC.Lon=latlon(2);
    % latlon=get_hdf_attribute(dset_id(i),'Top Right Geodetic Coordinates');
    % output_meta.GeoData.ImageCorners.ICP.LRFC.Lat=latlon(1);
    % output_meta.GeoData.ImageCorners.ICP.LRFC.Lon=latlon(2);

    %% Grid
    output_meta.Grid.Row.SS=get_hdf_attribute(dset_id(i),'Column Spacing');
    % output_meta.Grid.Row.SS=get_hdf_attribute(dset_id(i),'Column Time Interval')*SPEED_OF_LIGHT/2; % Exactly equivalent to above
    % Col.SS is derived after DRateSFPoly below, rather than used from this
    % given field, so that SICD metadata can be internally consistent:
    % output_meta.Grid.Col.SS=get_hdf_attribute(dset_id(i),'Line Spacing');
    output_meta.Grid.Row.ImpRespBW=2*get_hdf_attribute(group_id(i),'Range Focusing Bandwidth')/SPEED_OF_LIGHT;
    output_meta.Grid.Row.DeltaK1=-output_meta.Grid.Row.ImpRespBW/2;
    output_meta.Grid.Row.DeltaK2=-output_meta.Grid.Row.DeltaK1;
    if strcmpi(output_meta.Grid.Row.WgtType.WindowName,'HAMMING') % The usual CSM weigting
        % Computation of broadening factor for uniform window is:
        % 2 * fzero(@(x) (sin(pi*x)/(pi*x)) - (1/sqrt(2)), .1)
        a = str2double(output_meta.Grid.Row.WgtType.Parameter.value); % Generalized Hamming window parameter
        row_broadening_factor = 2*fzero(@(x) ...
            a*(sin(pi*x)/(pi*x)) + ((1-a)*(sin(pi*(x-1))/(pi*(x-1)))/2) + ...
            ((1-a)*(sin(pi*(x+1))/(pi*(x+1)))/2) - a/sqrt(2), .1);
        output_meta.Grid.Row.ImpRespWid = row_broadening_factor/output_meta.Grid.Row.ImpRespBW;
    end

    %% Timeline
    prf=get_hdf_attribute(group_id(i),'PRF');
    output_meta.Timeline.IPP.Set.IPPEnd=uint32(floor(prf*band_independent_meta.Timeline.CollectDuration));
    output_meta.Timeline.IPP.Set.IPPPoly=[0; prf];
    output_meta.Timeline.IPP.Set.TEnd=output_meta.Timeline.CollectDuration;

    %% RadarCollection
    % output_meta.RadarCollection.RefFreqIndex=uint32(0); % Absence of this field means all frequencies are true values
    chirp_length=get_hdf_attribute(group_id(i),'Range Chirp Length');
    chirp_rate=abs(get_hdf_attribute(group_id(i),'Range Chirp Rate'));
    bw=chirp_length*chirp_rate;
    output_meta.RadarCollection.TxFrequency.Min=fc-(bw/2);
    output_meta.RadarCollection.TxFrequency.Max=fc+(bw/2);
    output_meta.RadarCollection.Waveform.WFParameters.TxPulseLength=chirp_length;
    output_meta.RadarCollection.Waveform.WFParameters.TxRFBandwidth=bw;
    output_meta.RadarCollection.Waveform.WFParameters.TxFreqStart=...
        output_meta.RadarCollection.TxFrequency.Min;
    output_meta.RadarCollection.Waveform.WFParameters.TxFMRate=chirp_rate;
    sample_rate=get_hdf_attribute(group_id(i),'Sampling Rate');
    try
        if isnan(get_hdf_attribute(group_id(i),'Reference Dechirping Time'))
            output_meta.RadarCollection.Waveform.WFParameters.RcvDemodType='CHIRP';
            output_meta.RadarCollection.Waveform.WFParameters.RcvFMRate=0;
        else
            output_meta.RadarCollection.Waveform.WFParameters.RcvDemodType='STRETCH';
            output_meta.RadarCollection.Waveform.WFParameters.RcvFMRate=chirp_rate;
        end
    catch
        output_meta.RadarCollection.Waveform.WFParameters.RcvDemodType='CHIRP';
        output_meta.RadarCollection.Waveform.WFParameters.RcvFMRate=0;
    end
    output_meta.RadarCollection.Waveform.WFParameters.RcvWindowLength=...
        double(get_hdf_attribute(group_id(i),'Echo Sampling Window Length'))/sample_rate;
    output_meta.RadarCollection.Waveform.WFParameters.ADCSampleRate=sample_rate;

    %% ImageFormation
    output_meta.ImageFormation.RcvChanProc.ChanIndex=i;
    output_meta.ImageFormation.TxFrequencyProc.MinProc=...
        output_meta.RadarCollection.TxFrequency.Min;
    output_meta.ImageFormation.TxFrequencyProc.MaxProc=...
        output_meta.RadarCollection.TxFrequency.Max;
    if strcmp(mission_id, 'CSG')
        band_pol = deblank(get_hdf_attribute(HDF5_fid,'Polarization')');
        output_meta.ImageFormation.TxRcvPolarizationProc = [band_pol(1) ':' band_pol(2)];
    else
        output_meta.ImageFormation.TxRcvPolarizationProc=...
            band_independent_meta.RadarCollection.RcvChannels.ChanParameters(i).TxRcvPolarization;
    end

    %% RMA
    t_rg_scp = t_rg_first + (ss_rg_s*double(output_meta.ImageData.SCPPixel.Row)); % Range time to SCP
    output_meta.RMA.INCA.R_CA_SCP = t_rg_scp*SPEED_OF_LIGHT/2;
    t_az_scp = t_az_first + ...
        (ss_az_s*double(output_meta.ImageData.SCPPixel.Col)); % Zero doppler time of SCP
    % Compute DRateSFPoly
    % For the purposes of the DRateSFPoly computation, we ignore any
    % changes in velocity or doppler rate over the azimuth dimension.
    pos_coefs = [P_x(:) P_y(:) P_z(:)];
    % Velocity is derivate of position.
    vel_coefs=pos_coefs(1:end-1,:).*repmat(((size(pos_coefs,1)-1):-1:1)',[1 3]);
    scp_ca_time = t_az_scp + ref_time_offset;
    vel_x = polyval(vel_coefs(:,1), scp_ca_time);
    vel_y = polyval(vel_coefs(:,2), scp_ca_time);
    vel_z = polyval(vel_coefs(:,3), scp_ca_time);
    vm_ca_sq = vel_x.^2 + vel_y.^2 + vel_z.^2; % Magnitude of the velocity squared
    r_ca = [output_meta.RMA.INCA.R_CA_SCP; 1]; % Polynomial representing range as a function of range distance from SCP
    if exist('ss_rg_s','var')
        dop_rate_poly_rg_shifted=polyshift(dop_rate_poly_rg, t_rg_scp-t_rg_ref);
        dop_rate_poly_rg_scaled=dop_rate_poly_rg_shifted.*...
            (ss_rg_s/output_meta.Grid.Row.SS).^(0:(length(dop_rate_poly_rg)-1)).';
    end
    output_meta.RMA.INCA.DRateSFPoly = - conv(dop_rate_poly_rg_scaled,r_ca) * ... % Multiplication of two polynomials is just a convolution of their coefficients
        SPEED_OF_LIGHT / (2 * fc * vm_ca_sq(1)); % Assumes a SGN of -1
    % Fields dependent on Doppler rate
    output_meta.Grid.Col.SS = sqrt(vm_ca_sq(1)) * abs(ss_az_s) * output_meta.RMA.INCA.DRateSFPoly(1,1);
    % 'Azimuth Focusing Bandwidth' sometimes represents something larger
    % than ImpRespBW for S2 CSK mode.  Not sure what "transition bandwidth"
    % means, with respect to SAR, but it appears to correlate with non-zero
    % Doppler frequency content in SLC.
    switch mission_id
        case 'CSK'
            dop_bw = get_hdf_attribute(group_id(i),'Azimuth Focusing Transition Bandwidth'); % Doppler frequency
        otherwise
            dop_bw = get_hdf_attribute(group_id(i),'Azimuth Focusing Bandwidth'); % Doppler frequency
    end
    output_meta.Grid.Col.ImpRespBW = ... % Convert to azimuth spatial bandwidth (cycles per meter)
        min(dop_bw*abs(ss_az_s),1)/output_meta.Grid.Col.SS; % Can't have more bandwidth in data than sample spacing
    if strcmpi(output_meta.Grid.Col.WgtType.WindowName,'HAMMING') % The usual CSM weigting
        a = str2double(output_meta.Grid.Col.WgtType.Parameter.value); % Generalized Hamming window parameter
        col_broadening_factor = 2*fzero(@(x) ...
            a*(sin(pi*x)/(pi*x)) + ((1-a)*(sin(pi*(x-1))/(pi*(x-1)))/2) + ...
            ((1-a)*(sin(pi*(x+1))/(pi*(x+1)))/2) - a/sqrt(2), .1);
        output_meta.Grid.Col.ImpRespWid = col_broadening_factor/output_meta.Grid.Col.ImpRespBW;
    end
    output_meta.RMA.INCA.TimeCAPoly = [scp_ca_time; ... % With respect to start of collect
        ss_az_s/output_meta.Grid.Col.SS]; % Convert zero doppler spacing from sec/pixels to sec/meters
    % Compute DopCentroidPoly/DeltaKCOAPoly
    if exist('ss_rg_s','var')
        output_meta.RMA.INCA.DopCentroidPoly=zeros(length(dop_poly_rg),length(dop_poly_az));
        % Compute doppler centroid value at SCP
        output_meta.RMA.INCA.DopCentroidPoly(1) = ...
            polyval(dop_poly_rg(end:-1:1), t_rg_scp-t_rg_ref) + ...
            polyval(dop_poly_az(end:-1:1), t_az_scp-t_az_ref) - ...
            mean([dop_poly_az(1) dop_poly_rg(1)]); % These should be identical
        % Shift 1D polynomials to account for SCP
        dop_poly_az_shifted=polyshift(dop_poly_az, t_az_scp-t_az_ref);
        dop_poly_rg_shifted=polyshift(dop_poly_rg, t_rg_scp-t_rg_ref);
        % Scale 1D polynomials to from Hz/s^n to Hz/m^n
        dop_poly_az_scaled=dop_poly_az_shifted.'.*...
            (ss_az_s/output_meta.Grid.Col.SS).^(0:(length(dop_poly_az)-1));
        dop_poly_rg_scaled=dop_poly_rg_shifted.*...
            (ss_rg_s/output_meta.Grid.Row.SS).^(0:(length(dop_poly_rg)-1)).';
        output_meta.RMA.INCA.DopCentroidPoly(2:end,1)=dop_poly_rg_scaled(2:end);
        output_meta.RMA.INCA.DopCentroidPoly(1,2:end)=dop_poly_az_scaled(2:end);
        output_meta.RMA.INCA.DopCentroidCOA=true;
        output_meta.Grid.Col.DeltaKCOAPoly=...
            output_meta.RMA.INCA.DopCentroidPoly*ss_az_s/output_meta.Grid.Col.SS;
    end
    % TimeCOAPoly
    % TimeCOAPoly=TimeCA+(DopCentroid/dop_rate)
    % Since we can't evaluate this equation analytically, we will evaluate
    % samples of it across our image and fit a 2D polynomial to it.
    POLY_ORDER = 2; % Order of polynomial which we want to compute
    grid_samples = POLY_ORDER + 1; % in each dimension
    coords_az_m = linspace(-double(output_meta.ImageData.SCPPixel.Col),...
        double(output_meta.ImageData.NumCols-output_meta.ImageData.SCPPixel.Col-1), grid_samples) * ...
        output_meta.Grid.Col.SS;
    coords_rg_m = linspace(-double(output_meta.ImageData.SCPPixel.Row),...
        double(output_meta.ImageData.NumRows-output_meta.ImageData.SCPPixel.Row-1), grid_samples) * ...
        output_meta.Grid.Row.SS;
    timeca_sampled = sicd_polyval2d(output_meta.RMA.INCA.TimeCAPoly(:).',coords_az_m,coords_rg_m);
    dopcentroid_sampled = sicd_polyval2d(output_meta.RMA.INCA.DopCentroidPoly,coords_az_m,coords_rg_m);
    doprate_sampled = sicd_polyval2d(dop_rate_poly_rg_scaled.',coords_az_m,coords_rg_m);
    timecoapoly_sampled = timeca_sampled+(dopcentroid_sampled./doprate_sampled);
    % Least squares fit for 2D polynomial
    % A*x = b
    [coords_az_m, coords_rg_m] = ndgrid(coords_az_m, coords_rg_m);
    a = zeros(grid_samples^2, (POLY_ORDER+1)^2);
    for k = 0:POLY_ORDER
        for j = 0:POLY_ORDER
            a(:,k*(POLY_ORDER+1)+j+1) = (coords_rg_m(:).^j).*(coords_az_m(:).^k);
        end
    end
    b_coa = zeros((POLY_ORDER+1)^2,1);
    for k=1:((POLY_ORDER+1)^2)
       b_coa(k)=sum(timecoapoly_sampled(:).*a(:,k)); % center of aperture
    end
    A=zeros((POLY_ORDER+1)^2);
    for k=1:((POLY_ORDER+1)^2)
        for j=1:((POLY_ORDER+1)^2)
            A(k,j)=sum(a(:,k).*a(:,j));
        end
    end
    old_warning_state=warning('off','MATLAB:nearlySingularMatrix');
    x=A\b_coa; % MATLAB often flags this as badly scaled, but results still appear valid
    warning(old_warning_state);
    output_meta.Grid.TimeCOAPoly=reshape(x, POLY_ORDER+1, POLY_ORDER+1);

    %% Radiometric
    if strcmp(mission_id,'CSK')  % Unsure how to handle for CSG/KMPS
        if ~strncmpi(get_hdf_attribute(HDF5_fid,'Range Spreading Loss Compensation Geometry')', 'NONE', 4)
            fact = get_hdf_attribute(HDF5_fid,'Reference Slant Range') ^ ...
                (2 * get_hdf_attribute(HDF5_fid,'Reference Slant Range Exponent'));
            % This code commented out below converts from beta_0 to sigma_0,
            % but we will just populate the more natural beta_0 directly, and
            % let derived_sicd_fields compute the rest of the rcs, sigma_0,
            % gamma_0 terms in a sensor independent way later.
            % if ~strncmpi(get_hdf_attribute(HDF5_fid,'Incidence Angle Compensation Geometry')', 'NONE', 4)
            %     fact = fact * sind(get_hdf_attribute(HDF5_fid,'Reference Incidence Angle'));
            % end
            if isequal(get_hdf_attribute(HDF5_fid,'Calibration Constant Compensation Flag'), 0)
                fact = fact * (1 / (get_hdf_attribute(HDF5_fid,'Rescaling Factor')^2));
                fact = fact / get_hdf_attribute(group_id(i),'Calibration Constant');
                output_meta.Radiometric.BetaZeroSFPoly = fact;
            end
        end
    end

    %% GeoData
    % Now that sensor model fields have been populated, we can populate
    % GeoData.SCP more precisely.
    ecf = point_image_to_ground([output_meta.ImageData.SCPPixel.Row;output_meta.ImageData.SCPPixel.Col],output_meta);
    output_meta.GeoData.SCP.ECF.X=ecf(1);
    output_meta.GeoData.SCP.ECF.Y=ecf(2);
    output_meta.GeoData.SCP.ECF.Z=ecf(3);
    llh=ecf_to_geodetic([output_meta.GeoData.SCP.ECF.X output_meta.GeoData.SCP.ECF.Y output_meta.GeoData.SCP.ECF.Z]);
    output_meta.GeoData.SCP.LLH.Lat=llh(1);
    output_meta.GeoData.SCP.LLH.Lon=llh(2);
    output_meta.GeoData.SCP.LLH.HAE=llh(3);

    %% SCPCOA
    output_meta = derived_sicd_fields(output_meta);

    grouped_meta{i}=output_meta;
end
output_meta=grouped_meta;

%% Close HDF ids
for i=1:numbands
    H5S.close(dspace_id(i));
    H5D.close(dset_id(i));
    H5G.close(group_id(i));
end

end

% Reads and HDF5 attribute
function value = get_hdf_attribute(hid_t,attribute_name)
attribute=H5A.open_name(hid_t,attribute_name);
value = H5A.read(attribute,H5A.get_type(attribute));
H5A.close(attribute);
end

% Reproduces hamming functionality from the MATLAB Signal Processing
% Toolbox, but allows for arbitrary coefficients of raised cosine.
function w = raised_cos_fun(n, coef)
    w = coef - (1-coef)*cos(2*pi*(0:ceil(n/2)-1)'/(n-1));
    if ~rem(n,2)
        w = [w; w(end:-1:1)];
    else
        w = [w; w(end-1:-1:1)];
    end
end

% MATLAB's datenum function won't handle precise times down under a
% millisecond, because 1) It won't accept a format with more than 3 .FFF in
% the string description of the date format and 2) the resulting serial
% date number is stored in days from 00-JAN-0000 and just doesn't have
% enough bits to handle fractional seconds to the level we want.  Here we
% handle the fractional seconds separately so we can read date with the
% precision we need.
function [datenum_s, datenum_frac] = datenum_w_frac(datestring)
    datenum_s = datenum(datestring,'yyyy-mm-dd HH:MM:SS');
    datenum_frac = str2double(regexp(datestring,'\.\d*','match'));
    if isnan(datenum_frac), datenum_frac = 0; end;
end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////


text:
function qc=rotqr2qc(qr)
%ROTQR2QC converts a matrix of real quaternion vectors into complex form
%
% Inputs:
%
%     QR(4m,n)   mxn matrix of real-valued quaternions
%
% Outputs:
%
%     QC(2m,n)   mxn matrix of complex-valued quaternions
%
% The real-valued quaternion [r a b c]' becomes [r+j*b  a+j*c]'

%
%      Copyright (C) Mike Brookes 2000-2006
%      Version: $Id: rotqr2qc.m,v 1.3 2007/11/23 18:47:46 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=size(qr,1);
i=(1:2:n)-mod(0:n/2-1,2);
qc=qr(i,:)+1i*qr(i+2,:);


text:
function pf = max(pf1,pf2)
% maximum of two pole figures or the maximum of a single polefigure
%
% Syntax:
%   m  = max(pf)
%   pf = max(pf1,pf2)
%   pf = max(pf1,x) = max(x,pf1)
%
% Input
%  pf       - @PoleFigure
%  pf1, pf2 - @PoleFigure
%  x        - double
%
% Output
%  m  - double
%  pf - @PoleFigure

if nargin == 1

  pf = cellfun(@(x) max(x(:)),pf1.allI);

elseif isa(pf1,'double')

  pf = pf2;
  pf.intensities = max(pf1,pf2.intensities);

elseif isa(pf2,'double')

  pf = pf1;
  pf.intensities = max(pf2,pf1.intensities);

else

  pf = pf1;
  pf.intensities = max(pf1.intensities,pf2.intensities);

end

author: mtex-toolbox

repo: mtex

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


text:
function [theta, J_history] = gradientDescentMulti(X, y, theta, alpha, num_iters)
    %% GRADIENTDESCENTMULTI Performs gradient descent to learn theta
    %   theta = GRADIENTDESCENTMULTI(x, y, theta, alpha, num_iters) updates theta by
    %   taking num_iters gradient steps with learning rate alpha
    [theta, J_history] = gradientDescent(X, y, theta, alpha, num_iters);
end

author: worldveil

repo: coursera-ml


text:
function OK = convertCobraLP2mps(LPProblem, name)
% Creates an MPS (Mathematical Programming System) format ascii file
% representing the Linear Programming problem given by LPProblem.
%
% USAGE:
%
%    OK = convertCobraLP2mps(LPProblem, name)
%
% INPUT:
%    LPproblem:    Structure containing the following fields describing the LP problem to be solved
%
%                    * .A - LHS matrix
%                    * .b - RHS vector
%                    * .c - Objective coeff vector
%                    * .lb - Lower bound vector
%                    * .ub - Upper bound vector
%                    * .osense - Objective sense (max=-1, min=+1)
%                    * .csense - Constraint senses, a string containting the constraint sense for
%                      each row in A ('E', equality, 'G' greater than, 'L' less than).
%
% OPTIONAL INPUT:
%    name:         string giving name of LP problem
%
% OUTPUT:
%    OK:           1 if saving is success, 0 otherwise
%
% .. Authors:
%       - Ronan M.T. Fleming: 7 Sept 09
%       - Bruno Luong 03 Sep 2009  Uses MPS format exporting tool
%       http://www.mathworks.com/matlabcentral/fileexchange/19618
%
% ..
%    The MPS (Mathematical Programming System) file format was introduced by
%    IBM in 1970s, but has also been accepted by most subsequent linear
%    programming codes. To learn about MPS format, please see:
%    http://lpsolve.sourceforge.net/5.5/mps-format.htm

if ~exist('name', 'var')
    name = 'CobraLPProblem';
end

if isfield(LPProblem, 'S') && ~isfield(LPProblem, 'A')
    LPProblem.A = LPProblem.S;
end

mlt = size(LPProblem.A, 1);
if ~isfield(LPProblem, 'csense')
    LPProblem.csense(1:mlt) = 'E';
end
if size(LPProblem.csense, 1) > size(LPProblem.csense, 2)
    LPProblem.csense = LPProblem.csense';
end

E = false(mlt, 1);
G = false(mlt, 1);
L = false(mlt, 1);
Eind = findstr('E', LPProblem.csense);
Gind = findstr('G', LPProblem.csense);
Lind = findstr('L', LPProblem.csense);
E(Eind) = 1;
G(Gind) = 1;
L(Lind) = 1;

Aeq = LPProblem.A(E, :);
beq = LPProblem.b(E, 1);

% need to change sign of A*x >= b constraints
A2 = LPProblem.A;
b2 = LPProblem.b;
A2(G) = -A2(G);
b2(G) = -b2(G);
A = A2(G | L, :);
b = b2(G | L, :);

% Assume max c'v s.t. S v = b if osense not provided
if ~isfield(LPProblem, 'osense')
    LPProblem.osense = -1;
end

cost = LPProblem.c * LPProblem.osense;

L = LPProblem.lb;
U = LPProblem.ub;

% Build ascii fixed-width MPS matrix string that contains linear
% programming (LP) problem:
%
% Minimizing (for x in R^n): f(x) = cost'*x, subject to
%       A*x <= b        (LE)
%       Aeq*x = beq     (EQ)
%       L <= x <= U     (BD).

[Contain] = BuildMPS(A, b, Aeq, beq, cost, L, U, upper(name));

% Save matrix sring Contain in file "filename"
% Return OK == 1 if saving is success
%        OK == 0 otherwise
filename = [name '.mps'];
OK = SaveMPS(filename, Contain);

author: opencobra

repo: cobratoolbox


text:
function state = psocalculatepenalties(state)
% Calculates the penalty to apply for out-of-bounds particles when using
% the "penalize" constraint handling method.

% nx1 vector of penalties to apply to each particle
g = state.ConstrViolations ;
k = zeros(size(g));
n = size(state.Score,1);
% penalties = zeros(size(g,1),1) ;

if sum(state.OutOfBounds)
    idx = logical(state.OutOfBounds) ;
    k(idx,:) = (mean(state.Score).*g(idx,:))./sum(mean(g(idx,:),1).^2) ;
%     penalties = sum(k(idx,:).*g(idx,:),2) ;
    penalties = sum(k(idx,:).*g(idx,:),2) + max(state.Score(setdiff(1:n,idx))) ;

    % Apply the penalties to each particle's score
    state.Score(idx) = state.Score(idx) + penalties ;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function M = stiefelcomplexfactory(n, p, k)
% Returns a manifold struct. to optimize over complex orthonormal matrices.
%
% function M = stiefelcomplexfactory(n, p)
% function M = stiefelcomplexfactory(n, p, k)
%
% The complex Stiefel manifold is the set of complex orthonormal nxp
% matrices. If k is larger than 1, this is the Cartesian product of the
% complex Stiefel manifold taken k times. The metric is such that the
% manifold is a Riemannian submanifold of C^nxp equipped with the usual
% real-trace inner product, that is, it is the usual metric for the complex
% plane identified with R^2.
%
% Points are represented as matrices X of size n x p x k (or n x p if k=1,
% which is the default) such that each complex n x p matrix is orthonormal,
% i.e., X'*X = eye(p) if k = 1, or X(:, :, i)' * X(:, :, i) = eye(p) for
% i = 1 : k if k > 1. Tangent vectors are represented as matrices the same
% size as points.
%
% By default, k = 1.
%
%
% Please cite the Manopt paper as well as either of these research papers
% pertaining to this specific geometry:
% @InProceedings{sato2013complex,
%   Title        = {A complex singular value decomposition algorithm based on the {R}iemannian {N}ewton method},
%   Author       = {Sato, H. and Iwai, T.},
%   Booktitle    = {Decision and Control ({CDC}), 2013 {IEEE} 52nd Annual Conference on},
%   Year         = {2013},
%   Organization = {IEEE},
%   Pages        = {2972--2978}
% }
% @InProceedings{sato2014Riemannian,
%   Title        = {{R}iemannian conjugate gradient method for complex singular value decomposition problem},
%   Author       = {Sato, H.},
%   Booktitle    = {Decision and Control ({CDC}), 2014 {IEEE} 53rd Annual Conference on},
%   Year         = {2014},
%   Organization = {IEEE},
%   Pages        = {5849--5854}
% }
%
%
% See also: stiefelfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Hiroyuki Sato, April 27, 2015.
% Contributors:
% Change log:

    if ~exist('k', 'var') || isempty(k)
        k = 1;
    end

    if k == 1
        M.name = @() sprintf('Complex Stiefel manifold St(%d, %d)', n, p);
    elseif k > 1
        M.name = @() sprintf('Product complex Stiefel manifold St(%d, %d)^%d', n, p, k);
    else
        error('k must be an integer no less than 1.');
    end

    M.dim = @() k*(2*n*p - p^2); %! k*(n*p - .5*p*(p+1)) -> k*(2*n*p - p^2)

    M.inner = @(x, d1, d2) real(d1(:)'*d2(:)); %! trace -> real-trace

    M.norm = @(x, d) norm(d(:));

    M.dist = @(x, y) error('stiefel.dist not implemented yet.');

    M.typicaldist = @() sqrt(p*k);

    M.proj = @projection;
    function Up = projection(X, U)

        XHU = multiprod(multihconj(X), U); %! XtU -> XHU, multitransp -> multihconj
        herXHU = multiherm(XHU); %! symXtU -> herXHU, multisym -> multiherm
        Up = U - multiprod(X, herXHU); %! symXtU -> herXHU

    end

    M.tangent = M.proj;

    % For Riemannian submanifolds, converting a Euclidean gradient into a
    % Riemannian gradient amounts to an orthogonal projection.
	M.egrad2rgrad = M.proj;

    M.ehess2rhess = @ehess2rhess;
    function rhess = ehess2rhess(X, egrad, ehess, H)
        XHG = multiprod(multihconj(X), egrad); %! XtG -> XHG, multitransp -> multihconj
        herXHG = multiherm(XHG); %! symXtG -> herXHG, multisym(XtG) -> multiherm(XHG)
        HherXHG = multiprod(H, herXHG); %! HsymXtG -> HherXHG, symXtG -> herXHG
        rhess = projection(X, ehess - HherXHG); %! HsymXtG -> HherXHG
    end

    M.retr = @retraction;
    function Y = retraction(X, U, t)
        if nargin < 3
            t = 1.0;
        end
        Y = X + t*U;
        for i = 1 : k
            [Q, R] = qr(Y(:, :, i), 0);
            % The instruction with R assures we are not flipping signs
            % of some columns, which should never happen in modern Matlab
            % versions but may be an issue with older versions.
            Y(:, :, i) = Q * diag(sign(sign(diag(R))+.5));
        end
    end

    M.exp = @exponential;
    function Y = exponential(X, U, t)
        if nargin == 2
            t = 1;
        end
        tU = t*U;
        Y = zeros(size(X));
        for i = 1 : k
            % From a formula by Ross Lippert, Example 5.4.2 in AMS08.
            Xi = X(:, :, i);
            Ui = tU(:, :, i);
            Y(:, :, i) = [Xi Ui] * ...
                         expm([Xi'*Ui , -Ui'*Ui ; eye(p) , Xi'*Ui]) * ...
                         [ expm(-Xi'*Ui) ; zeros(p) ];
        end

    end

    M.hash = @(X) ['z' hashmd5([real(X(:)) ; imag(X(:))])]; %! X(:) -> [real(X(:)) ; imag(X(:))]

    M.rand = @random;
    function X = random()
        X = zeros(n, p, k);
        for i = 1 : k
            [Q, unused] = qr(randn(n, p) + 1i*randn(n,p), 0); %#ok<NASGU> %! Complex version
            X(:, :, i) = Q;
        end
    end

    M.randvec = @randomvec;
    function U = randomvec(X)
        U = projection(X, randn(n, p, k) + 1i*randn(n, p, k)); %! Complex version
        U = U / norm(U(:));
    end

    M.lincomb = @matrixlincomb;

    M.zerovec = @(x) zeros(n, p, k);

    M.transp = @(x1, x2, d) projection(x2, d);

    M.vec = @(x, u_mat) [real(u_mat(:)) ; imag(u_mat(:))];
    M.mat = @(x, u_vec) reshape(u_vec(1:(n*p*k)) + 1i*u_vec((n*p*k+1):end), [n, p, k]);
    M.vecmatareisometries = @() true; % TODO : to check.

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


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

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

% Check if any coefficients are NaN:
out = any(isnan(f.coeffs(:)));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function V = LP10(K, P, v, LPproblem, epsilon)
% Finds a flux vector that maintains the activity of any active irreversible
% core reaction (K) yet minimises the activity of any non-core reaction (P).
% Implementation of LP-10 for input sets K, P (see FASTCORE paper)
%
% USAGE:
%
%    V = LP10(K, P, v, epsilon)
%

% .. Authors: -  Nikos Vlassis, Maria Pires Pacheco, Thomas Sauter, 2013
%             LCSB / LSRU, University of Luxembourg
%   2019/04/08: Agnieszka Wegrzyn - updated the function to work with models with coupling constraints
% 2020 Ronan Fleming - adaptive scaling prototype

V = [];
if isempty(P) || isempty(K)
    return;
end
[m2,n2] = size(LPproblem.A);

%force active irreversible core reactions to be active
%important to set lower bound to min(v(K))*10 not min(v(K)) as the latter is too close to the threshold
minIrrevCoreSetBound = min(v(K))*10;
%disp(minIrrevCoreSetBound);

np = numel(P);
nk = numel(K);
%[m,n] = size(model.S);
[m2,n2] = size(LPproblem.A);

% objective
c = [zeros(n2,1); ones(np,1)];

% S*v = b and C *v<= d if present
Aeq = [LPproblem.A, sparse(m2,np)]; %changed the size of sparse() to match the size of LPproblem.A
beq = LPproblem.b;

% inequalities
Ip = sparse(np,n2); Ip(sub2ind(size(Ip),(1:np)',P(:))) = 1;
Ik = sparse(nk,n2); Ik(sub2ind(size(Ik),(1:nk)',K(:))) = 1;
Aineq = sparse(...
    [[Ip, -speye(np)]; ...  %  v - z <= 0                         v \in P
    [-Ip, -speye(np)]; ...  % -v - z <= 0                         v \in P
    [-Ik, sparse(nk,np)]]); % -v     <= -minIrrevCoreSetBound     v \in K
bineq = [zeros(2*np,1); -ones(nk,1)*minIrrevCoreSetBound];

% bounds
vlb = LPproblem.lb;
%vlb(K)=minIrrevCoreSetBound;
lb = [vlb; zeros(np,1)];
ub = [LPproblem.ub; max(abs(LPproblem.ub(P)),abs(LPproblem.lb(P)))];

% Set up LP problem
LP10problem.A=[Aeq;Aineq];
LP10problem.b=[beq;bineq];
LP10problem.lb=lb;
LP10problem.ub=ub;
LP10problem.c=c;
LP10problem.osense=1;%minimise
LP10problem.csense = [LPproblem.csense; repmat('L',2*np + nk,1)];

solution = solveCobraLP(LP10problem);
if solution.stat~=1
    fprintf('%s\n','LP10: Adaptive scaling factor did not work, switching to fixed scaling factor of 1e4')

    scalingfactor = 1e4;

    np = numel(P);
    nk = numel(K);

    % objective
    f = [zeros(n2,1); ones(np,1)];

    % S*v = b and C *v<= d if present
    Aeq = [LPproblem.A, sparse(m2,np)]; %changed the size of sparse() to match the size of LPproblem.A
    beq = LPproblem.b;

    % inequalities
    Ip = sparse(np,n2); Ip(sub2ind(size(Ip),(1:np)',P(:))) = 1;
    Ik = sparse(nk,n2); Ik(sub2ind(size(Ik),(1:nk)',K(:))) = 1;
    Aineq = sparse([[Ip, -speye(np)]; ...
        [-Ip, -speye(np)]; ...
        [-Ik, sparse(nk,np)]]);
    bineq = [zeros(2*np,1); -ones(nk,1)*epsilon];

    % bounds
    lb = [LPproblem.lb; zeros(np,1)];
    ub = [LPproblem.ub; max(abs(LPproblem.ub(P)),abs(LPproblem.lb(P)))] ;

    % Set up LP problem
    LP10problem.A=[Aeq;Aineq];
    %scaling factor must be multiplied by vector on rhs of A*x<= b
    LP10problem.b=[beq;bineq]*scalingfactor;
    %scaling factor must be multiplied by lb and ub of lb <= x <= ub
    LP10problem.lb=lb*scalingfactor;
    LP10problem.ub=ub*scalingfactor;
    LP10problem.c=f;
    LP10problem.osense=1;%minimise
    LP10problem.csense = [LPproblem.csense; repmat('L',2*np + nk,1)];

    solution = solveCobraLP(LP10problem);
    if solution.stat~=1
        disp(solution)
        warning(['solution.stat~=1. Scaling factor in LP10 is ' num2str(scalingfactor) ', which may cause numerical instability for solveCobraLP.'])
    end
end

if solution.stat~=1
    fprintf('\n%s%s\n',num2str(solution.stat),' = sol.stat')
    fprintf('%s%s\n',num2str(solution.origStat),' = sol.origStat')
    warning('LP solution may not be optimal')
end

x=solution.full;

if ~isempty(x)
    V = x(1:n2);
else
    V=ones(n2,1)*NaN;
end


text:
function [dat] = read_combined_ds(dirname, hdr, begsample, endsample, chanindx)

% READ_COMBINED_DS
%
% Use as
%   hdr = read_combined_ds(dirname)
%   dat = read_combined_ds(dirname, hdr, begsample, endsample, chanindx)

% Copyright (C) 2008, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: read_combined_ds.m 2892 2011-02-16 13:34:39Z crimic $

needhdr = nargin==1 || isempty(hdr);
needdat = nargin>1;

supported = {
  'plexon_nex'
  'neuralynx_bin'
  'neuralynx_ncs'
  'fcdc_matbin'
  };

if needhdr

  ls = dir(dirname);                % get the list of filenames
  ls = ls(~cell2mat({ls.isdir}));   % remove the subdirs

  nfiles = length(ls);
  fname = cell(nfiles,1);
  ftype = cell(nfiles,1);
  sel   = false(nfiles,1);

  for i=1:nfiles
    fname{i} = fullfile(dirname, ls(i).name);
    ftype{i} = ft_filetype(fname{i});
    sel(i)   = any(strcmp(ftype{i}, supported));
    [p, f, x] = fileparts(fname{i});
    if filetype_check_extension(fname{i}, '.mat')
      % select only the *.bin and not the *.mat of each pair
      sel(i) = false;
    end
  end

  if ~any(sel)
    error('no supported files were found');
  end

  fname = fname(sel);
  ftype = ftype(sel);
  nfiles = length(fname);

  clear filehdr
  for i=1:nfiles
    filehdr(i) = ft_read_header(fname{i}, 'headerformat', ftype{i});
  end

  if any([filehdr.nChans]~=1)
    error('more than one channel per file not supported');
  else
    hdr.nChans = sum([filehdr.nChans]);
  end

  if length(unique([filehdr.label]))~=nfiles
    error('not all channels have a unique name');
  else
    hdr.label = [filehdr.label]';
  end

  if any(diff([filehdr.Fs]))
    error('different sampling frequenties per file not supported');
  else
    hdr.Fs = filehdr(1).Fs;
  end

  if any(diff([filehdr.nSamples]))
    error('different nSamples per file not supported');
  else
    hdr.nSamples = filehdr(1).nSamples;
  end

  if any(diff([filehdr.nSamplesPre]))
    error('different nSamplesPre per file not supported');
  else
    hdr.nSamplesPre = filehdr(1).nSamplesPre;
  end

  if any(diff([filehdr.nTrials]))
    error('different nTrials per file not supported');
  else
    hdr.nTrials = filehdr(1).nTrials;
  end

  if isfield(filehdr, 'TimeStampPerSample') && any(diff([filehdr.TimeStampPerSample]))
    error('different TimeStampPerSample per file not supported');
  else
    hdr.TimeStampPerSample = filehdr(1).TimeStampPerSample;
  end

  if isfield(filehdr, 'FirstTimeStamp') && any(diff([filehdr.FirstTimeStamp]))
    error('different FirstTimeStamp per file not supported');
  else
    hdr.FirstTimeStamp = filehdr(1).FirstTimeStamp;
  end

  % remember the original header details
  hdr.orig.header = filehdr;
  hdr.orig.fname  = fname;
  hdr.orig.ftype  = ftype;
end

if needdat
  filehdr = hdr.orig.header;
  fname   = hdr.orig.fname;
  ftype   = hdr.orig.ftype;

  if nargin<3
    begsample=1;
  end
  if nargin<4
    endsample = hdr.nSamples*hdr.nTrials;
  end
  if nargin<5
    chanindx = 1:hdr.nChans;
  end

  nsamples  = endsample-begsample+1;
  nchans    = length(chanindx);
  dat       = zeros(nchans, nsamples);

  for i=1:nchans
    thischan = chanindx(i);
    tmp = read_data(fname{thischan}, 'header', filehdr(thischan), 'dataformat', ftype{thischan}, 'begsample', begsample, 'endsample', endsample);
    dat(i,:) = tmp;
  end
end

if ~needdat
  % return the header
  dat = hdr;
end

author: PatternRecognition

repo: OpenBMI


text:
function [bus, gen, branch] = pfsoln(baseMVA, bus0, gen0, branch0, Ybus, Yf, Yt, V, ref, pv, pq, mpopt)
%PFSOLN  Updates bus, gen, branch data structures to match power flow soln.
%   [BUS, GEN, BRANCH] = PFSOLN(BASEMVA, BUS0, GEN0, BRANCH0, ...
%                                   YBUS, YF, YT, V, REF, PV, PQ, MPOPT)

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

%% define named indices into bus, gen, branch matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;

%% default options
if nargin < 12
    mpopt = mpoption();
end

%% initialize return values
bus     = bus0;
gen     = gen0;
branch  = branch0;

%%----- update bus voltages -----
bus(:, VM) = abs(V);
bus(:, VA) = angle(V) * 180 / pi;

%%----- update Qg for gens at PV/slack buses and Pg for slack bus(es) -----
%% generator info
on = find(gen(:, GEN_STATUS) > 0 & ...  %% which generators are on?
        bus(gen(:, GEN_BUS), BUS_TYPE) ~= PQ);  %% ... and not at PQ buses
off = find(gen(:, GEN_STATUS) <= 0);    %% which generators are off?
gbus = gen(on, GEN_BUS);                %% what buses are they at?

%% compute total injected bus powers
Sbus = V(gbus) .* conj(Ybus(gbus, :) * V);

%% update Qg for generators at PV/slack buses
gen(off, QG) = zeros(length(off), 1);   %% zero out off-line Qg
%% don't touch the ones at PQ buses
[Pd_gbus, Qd_gbus] = total_load(bus(gbus, :), [], 'bus', [], mpopt);
gen(on, QG) = imag(Sbus) * baseMVA + Qd_gbus;   %% inj Q + local Qd
%% ... at this point any buses with more than one generator will have
%% the total Q dispatch for the bus assigned to each generator. This
%% must be split between them. We do it first equally, then in proportion
%% to the reactive range of the generator.

if length(on) > 1
    %% build connection matrix, element i, j is 1 if gen on(i) at bus j is ON
    nb = size(bus, 1);
    ngon = size(on, 1);
    Cg = sparse((1:ngon)', gbus, ones(ngon, 1), ngon, nb);

    %% divide Qg by number of generators at the bus to distribute equally
    ngg = Cg * sum(Cg)';    %% ngon x 1, number of gens at this gen's bus
    gen(on, QG) = gen(on, QG) ./ ngg;

    %% set finite proxy M for infinite limits (for ~ proportional splitting)
    %% equal to sum over all gens at bus of abs(Qg) plus any finite Q limits
    Qmin = gen(on, QMIN);
    Qmax = gen(on, QMAX);
    M = abs(gen(on, QG));
    M(~isinf(Qmax)) = M(~isinf(Qmax)) + abs(Qmax(~isinf(Qmax)));
    M(~isinf(Qmin)) = M(~isinf(Qmin)) + abs(Qmin(~isinf(Qmin)));
    M = Cg * Cg' * M;   %% each gen gets sum over all gens at same bus
    %% replace +/- Inf limits with proxy +/- M
    Qmin(Qmin ==  Inf) =  M(Qmin ==  Inf);
    Qmin(Qmin == -Inf) = -M(Qmin == -Inf);
    Qmax(Qmax ==  Inf) =  M(Qmax ==  Inf);
    Qmax(Qmax == -Inf) = -M(Qmax == -Inf);

    %% divide proportionally
    Cmin = sparse((1:ngon)', gbus, Qmin, ngon, nb);
    Cmax = sparse((1:ngon)', gbus, Qmax, ngon, nb);
    Qg_tot = Cg' * gen(on, QG);     %% nb x 1 vector of total Qg at each bus
    Qg_min = sum(Cmin)';            %% nb x 1 vector of min total Qg at each bus
    Qg_max = sum(Cmax)';            %% nb x 1 vector of max total Qg at each bus
    gen(on, QG) = Qmin + ...
        (Cg * ((Qg_tot - Qg_min)./(Qg_max - Qg_min + eps))) .* ...
            (Qmax - Qmin);          %%                ^ avoid div by 0

    %% fix gens at buses with Qg range = 0 (use equal violation for all)
    ig = find(abs(Cg * (Qg_min - Qg_max)) < 10*eps);  %% gens at buses with Qg range = 0
    if ~isempty(ig)
        ib = find(sum(Cg(ig,:), 1)');   %% buses with Qg range = 0
        %% total mismatch @ bus div by number of gens
        mis = sparse(ib, 1, (Qg_tot(ib) - Qg_min(ib)) ./ sum(Cg(:, ib)', 2), nb, 1);
        gen(on(ig), QG) = Qmin(ig) + Cg(ig, :) * mis;
    end
end                                             %% (terms are mult by 0 anyway)

%% update Pg for slack gen(s)
for k = 1:length(ref)
    refgen = find(gbus == ref(k));              %% which is(are) the reference gen(s)?
    Pd_refk = total_load(bus(ref(k), :), [], 'bus', [], mpopt);
    gen(on(refgen(1)), PG) = real(Sbus(refgen(1))) * baseMVA + Pd_refk; %% inj P + local Pd
    if length(refgen) > 1       %% more than one generator at this ref bus
        %% subtract off what is generated by other gens at this bus
        gen(on(refgen(1)), PG) = gen(on(refgen(1)), PG) ...
                                - sum(gen(on(refgen(2:length(refgen))), PG));
    end
end

%%----- update/compute branch power flows -----
out = find(branch(:, BR_STATUS) == 0);      %% out-of-service branches
br = find(branch(:, BR_STATUS));            %% in-service branches
Sf = V(branch(br, F_BUS)) .* conj(Yf(br, :) * V) * baseMVA; %% complex power at "from" bus
St = V(branch(br, T_BUS)) .* conj(Yt(br, :) * V) * baseMVA; %% complex power injected at "to" bus
branch(br, [PF, QF, PT, QT]) = [real(Sf) imag(Sf) real(St) imag(St)];
branch(out, [PF, QF, PT, QT]) = zeros(length(out), 4);

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%% Ex. 7 Elementary functions with a vectorial variable

a = [2 3 5];
b = sin(a)
% Output:
%        b = 0.9092 0.1411 -0.9589

%Remark: The content of b is [sin(2) sin(3) sin(5)].

author: TheAlgorithms

repo: MATLAB-Octave


text:
classdef VBoxFlex < uix.VBoxFlex
    %uiextras.VBoxFlex  A dynamically resizable vertical layout
    %
    %   obj = uiextras.VBoxFlex() creates a new dynamically resizable
    %   vertical box layout with all parameters set to defaults. The output
    %   is a new layout object that can be used as the parent for other
    %   user-interface components.
    %
    %   obj = uiextras.VBoxFlex(param,value,...) also sets one or more
    %   parameter values.
    %
    %   See the documentation for more detail and the list of properties.
    %
    %   Examples:
    %   >> f = figure( 'Name', 'uiextras.VBoxFlex example' );
    %   >> b = uiextras.VBoxFlex( 'Parent', f );
    %   >> uicontrol( 'Parent', b, 'Background', 'r' )
    %   >> uicontrol( 'Parent', b, 'Background', 'b' )
    %   >> uicontrol( 'Parent', b, 'Background', 'g' )
    %   >> uicontrol( 'Parent', b, 'Background', 'y' )
    %   >> set( b, 'Sizes', [-1 100 -2 -1], 'Spacing', 5 );
    %
    %   See also: uiextras.HBoxFlex
    %             uiextras.VBox
    %             uiextras.Grid

    %  Copyright 2009-2014 The MathWorks, Inc.
    %  $Revision: 979 $ $Date: 2014-09-28 14:26:12 -0400 (Sun, 28 Sep 2014) $

    properties( Hidden, Access = public, Dependent )
        Enable % deprecated
        Sizes
        MinimumSizes
        ShowMarkings
    end

    methods

        function obj = VBoxFlex( varargin )

            % Call uix constructor
            obj@uix.VBoxFlex( varargin{:} )

            % Auto-parent
            if ~ismember( 'Parent', varargin(1:2:end) )
                obj.Parent = gcf();
            end

        end % constructor

    end % structors

    methods

        function value = get.Enable( ~ )

            % Warn
            % warning( 'uiextras:Deprecated', ...
            %     'Property ''Enable'' will be removed in a future release.' )

            % Return
            value = 'on';

        end % get.Enable

        function set.Enable( ~, value )

            % Check
            assert( ischar( value ) && any( strcmp( value, {'on','off'} ) ), ...
                'uiextras:InvalidPropertyValue', ...
                'Property ''Enable'' must be ''on'' or ''off''.' )

            % Warn
            % warning( 'uiextras:Deprecated', ...
            %     'Property ''Enable'' will be removed in a future release.' )

        end % set.Enable

        function value = get.Sizes( obj )

            % Get
            value = transpose( obj.Heights );

        end % get.Sizes

        function set.Sizes( obj, value )

            % Set
            obj.Heights = value;

        end % set.Sizes

        function value = get.MinimumSizes( obj )

            % Get
            value = transpose( obj.MinimumHeights );

        end % get.MinimumSizes

        function set.MinimumSizes( obj, value )

            % Get
            obj.MinimumHeights = value;

        end % set.MinimumSizes

        function value = get.ShowMarkings( obj )

            % Get
            value = obj.DividerMarkings;

        end % get.ShowMarkings

        function set.ShowMarkings( obj, value )

            % Set
            obj.DividerMarkings = value;

        end % set.ShowMarkings

    end % accessors

end % classdef

author: goGPS-Project