text:
function plotExperimentalvsPredictedExchange(measuredFluxes, predictedFluxes, otherFluxes, extraFluxes, param)
% plot a comparison of experimental and predicted exchange reaction rates
%
% INPUT
% measuredFluxes:   table with the fluxes obtained from exometabolomics experiments with the following variables
%           * .rxns:    k x 1 cell array of reaction identifiers
%           * .mean:    k x 1 double of measured mean flux
%           * .SD:      k x 1 double standard deviation of the measured flux
%           * .labels:  k x 1 cell array of labels to display on y axis
%           * .Properties.Description: string describing the data, used in plot legend
%
% predictedFluxes:   table with the following variables
%           * .rxns:    k x 1 cell array of reaction identifiers
%           * .v:       k x 1 double of predicted mean flux
%           * .lb:      k x 1 double lower bound on reaction flux
%           * .ub:      k x 1 double upper bound on reaction flux
%           * .labels:  k x 1 cell array of labels to display on y axis
%           * .Properties.Description: string describing the data, used in plot legend
%
% OPTIONAL INPUT
% otherFluxes:  table with the following variables
%           * .rxns:    k x 1 cell array of reaction identifiers
%           * .Properties.Description: sring describing the data, used in plot legend
%           * .labels:  k x 1 cell array of labels to display on y axis
%   and
%           * .mean:    k x 1 double of measured mean flux
%           * .SD:      k x 1 double standard deviation of the measured flux
%
%   or
%           * .v:       k x 1 double of predicted mean flux
%           * .lb:      k x 1 double lower bound on reaction flux
%           * .ub:      k x 1 double upper bound on reaction flux
%
%
% param: paramters structure with the following fields
%           * .measuredFluxes  {(1),0} set to zero to not display experimental fluxes
%
%           * .labelType determines which y axes lables to display
%                       'metabolitePlatform' = platform left and metabolite name right
%                       'metabolite'         = metabolite only right  (default)
%                       'platform'           = platform only left
%           * .expOrder {(1),0}
%                       1  = order reactions by value of experimental flux
%                       0  = order reactions by value of predicted flux (predictedFluxes) for each reaction corresponding to an experimental flux
%           * .saveFigures = 1 saves figure as .fig and .png to current directory
%           * .plotTitle: title of the plot

if exist('measuredFluxes','var')
    expDescription = measuredFluxes.Properties.Description;
else
    error('predictedFluxes missing: must provide a predicted solution');
end

if exist('predictedFluxes','var')
    predDescription = predictedFluxes.Properties.Description;
else
    error('predictedFluxes missing: must provide a predicted solution');
end
if exist('otherFluxes','var') && ~isempty(otherFluxes)
    otherDescription = otherFluxes.Properties.Description;
else
    otherFluxes = [];
    otherDescription = '';
end
if exist('extraFluxes','var') && ~isempty(extraFluxes)
    extraDescription = extraFluxes.Properties.Description;
else
    extraFluxes = [];
    extraDescription = '';
end
if ~exist('param','var')
    param = struct();
end
if ~isfield(param,'measuredFluxes')
    param.measuredFluxes = 1;
end
if ~isfield(param,'labelType')
    param.labelType = 'metabolite';
end
if ~isfield(param,'saveFigures')
    param.saveFigures = 0;
end
if ~isfield(param,'expOrder')
    param.expOrder = 1;
end
if ~isfield(param,'plotTitle')
    param.plotTitle = [];
end
if ~isfield(param,'measuredBounds')
    param.measuredBounds = 1;
end
if ~isfield(param,'predictedBounds')
    param.predictedBounds = 1;
end
if ~isfield(param,'otherBounds')
    param.otherBounds = 1;
end
if ~isfield(param,'extraBounds')
    param.extraBounds = 1;
end
if ~isvar(measuredFluxes,'labels')
    if isvar(measuredFluxes,'metNames')
        measuredFluxes.labels = measuredFluxes.metNames;
    end
end
if isvar(measuredFluxes,'rxns')
    measuredFluxes.rxns = measuredFluxes.rxns;
end
if isvar(predictedFluxes,'rxns')
    predictedFluxes.rxns = predictedFluxes.rxns;
end
if ~isempty(otherFluxes)
    if isvar(otherFluxes,'rxns')
        otherFluxes.rxns = otherFluxes.rxns;
    end
    if ~isvar(otherFluxes,'rxns')
        error('otherFluxes.rxns missing - needed to map to other variables')
    end
end
if ~isempty(extraFluxes)
    if isvar(extraFluxes,'rxns')
        extraFluxes.rxns = extraFluxes.rxns;
    end
    if ~isvar(extraFluxes,'rxns')
        error('extraFluxes.rxns missing - needed to map to other variables')
    end
end
%%%%%
if param.expOrder
    %eliminate NaN experimental data
    measuredFluxes = measuredFluxes(~isnan(measuredFluxes.mean),:);

    %sort by the normal cumulative distribution function that the flux is an uptake (1) rather than a secretion (0)
    zerobool = (measuredFluxes.mean-measuredFluxes.SD)<=0 & (measuredFluxes.mean+measuredFluxes.SD)>=0;
    posbool = sign(measuredFluxes.mean)==1 & ~zerobool;
    negbool = sign(measuredFluxes.mean)==-1 & ~zerobool;

    measuredFluxes.uptakeProbability = normcdf(0,measuredFluxes.mean,measuredFluxes.SD,'upper');

    measuredFluxes.toRankOrder = zeros(length(measuredFluxes.mean),1);
    measuredFluxes.toRankOrder(zerobool) = measuredFluxes.uptakeProbability(zerobool);
    measuredFluxes.toRankOrder(posbool) =  measuredFluxes.mean(posbool)+1;
    measuredFluxes.toRankOrder(negbool) = measuredFluxes.mean(negbool);
    measuredFluxes.zerobool=zerobool;
    measuredFluxes.posbool=posbool;
    measuredFluxes.negbool=negbool;

    if 1
        [sortedUptakeProbability,xi] = sort(measuredFluxes.toRankOrder);
    else
        [sortedUptakeProbability,xi] = sort(measuredFluxes.uptakeProbability,'descend');
    end

    measuredFluxes = measuredFluxes(xi,:);
    zerobool = measuredFluxes.zerobool;
    posbool = measuredFluxes.posbool;
    negbool = measuredFluxes.negbool;

    %place the predicted exchanges in the same order as the measurements
    predictedFluxes = mapAontoB(predictedFluxes.rxns,measuredFluxes.rxns, predictedFluxes);
else
    for i=1:length(predictedFluxes.labels)
        labels2 = strrep(predictedFluxes.labels{i},' for ', ' of ');
        predictedFluxes.labels{i} = [predictedFluxes.labels{i} ', ' labels2];
    end

    %order the predicted exchanges
    predictedFluxes.zerobool = (predictedFluxes.v-1e-5)<=0 & (predictedFluxes.v+1e-5)>=0;
    predictedFluxes.posbool = sign(predictedFluxes.v)==1 & ~predictedFluxes.zerobool;
    predictedFluxes.negbool = sign(predictedFluxes.v)==-1 & ~predictedFluxes.zerobool;
    predictedFluxes.uptakeProbability = normcdf(0,predictedFluxes.v,abs(predictedFluxes.v)/10);

    toRankOrder = zeros(length(predictedFluxes.v),1);
    toRankOrder(predictedFluxes.zerobool) = predictedFluxes.uptakeProbability(predictedFluxes.zerobool);
    toRankOrder(predictedFluxes.posbool) =  predictedFluxes.v(predictedFluxes.posbool)+1;
    toRankOrder(predictedFluxes.negbool)=predictedFluxes.v(predictedFluxes.negbool);

    if 1
        [sortedUptakeProbability,xi] = sort(toRankOrder);
    else
        [sortedUptakeProbability,xi] = sort(predictedFluxes.uptakeProbability,'descend');
    end
    predictedFluxes = predictedFluxes(xi,:);
    zerobool = predictedFluxes.zerobool;
    posbool = predictedFluxes.posbool;
    negbool = predictedFluxes.negbool;

    %place the measurements in the same order as the predictions
    measuredFluxes = mapAontoB(measuredFluxes.rxns,predictedFluxes.rxns,measuredFluxes);
end

%log modulus transformation
x = logmod(measuredFluxes.mean,10);
xl = logmod(measuredFluxes.mean - measuredFluxes.SD,10);
xl = x - xl;
xu = logmod(measuredFluxes.mean + measuredFluxes.SD,10);
xu = xu - x;

%use log modulus transformation
v = logmod(predictedFluxes.v,10);
if isvar(predictedFluxes,'lb') && isvar(predictedFluxes,'ub')
    vl = logmod(predictedFluxes.lb,10);
    vl = v - vl;
    vu = logmod(predictedFluxes.ub,10);
    vu = vu - v;
else
    xl = logmod(measuredFluxes.mean - measuredFluxes.SD,10);
    xl = x - xl;
    xu = logmod(measuredFluxes.mean + measuredFluxes.SD,10);
    xu = xu - x;
end

if isempty(otherFluxes)
    if param.measuredFluxes
        min_v = min([v;x-xl]);
        max_v = max([v;x+xu]);
    else
        min_v = min(v);
        max_v = max(v);
    end
    bool=~isnan(v);
else
    if param.expOrder
        %place the other fluxes in the same order as the measurements
        otherFluxes = mapAontoB(otherFluxes.rxns,measuredFluxes.rxns, otherFluxes);
    else
        %place the other fluxes in the same order as the predictions
        otherFluxes = mapAontoB(otherFluxes.rxns,predictedFluxes.rxns, otherFluxes);
    end

    if isvar(otherFluxes,'v') && isvar(otherFluxes,'lb') && isvar(otherFluxes,'ub')
        %use log modulus transformation
        v2 = logmod(otherFluxes.v,10);
        v2l = logmod(otherFluxes.lb,10);
        v2l = v2 - v2l;
        v2u = logmod(otherFluxes.ub,10);
        v2u = v2u - v2;
    elseif isvar(otherFluxes,'mean') && isvar(otherFluxes,'SD')
        %use log modulus transformation
        v2 = logmod(otherFluxes.mean,10);
        v2l = logmod(otherFluxes.mean - otherFluxes.SD,10);
        v2l = v2 - v2l;
        v2u = logmod(otherFluxes.mean + otherFluxes.SD,10);
        v2u = v2u - v2;
    else
        error('Unrecognised input: otherFluxes')
    end

    %limits
    if param.measuredFluxes
        min_v = min([v;v2;x-xl]);
        max_v = max([v;v2;x+xu]);
    else
        min_v = min([v;v2]);
        max_v = max([v;v2]);
    end
    bool=~isnan(v) & ~isnan(v2);
end

if isempty(extraFluxes)
    if isempty(otherFluxes)
        if param.measuredFluxes
            min_v = min([v;x-xl]);
            max_v = max([v;x+xu]);
        else
            min_v = min(v);
            max_v = max(v);
        end
        bool=~isnan(v);
    else
        if param.measuredFluxes
            min_v = min([v;v2;x-xl]);
            max_v = max([v;v2;x+xu]);
        else
            min_v = min([v;v2]);
            max_v = max([v;v2]);
        end
        bool=~isnan(v) & ~isnan(v2);
    end
else
    if param.expOrder
        %place the other fluxes in the same order as the measurements
        extraFluxes = mapAontoB(extraFluxes.rxns,measuredFluxes.rxns, extraFluxes);
    else
        %place the other fluxes in the same order as the predictions
        extraFluxes = mapAontoB(extraFluxes.rxns,predictedFluxes.rxns, extraFluxes);
    end

    if isvar(extraFluxes,'v') && isvar(extraFluxes,'lb') && isvar(extraFluxes,'ub')
        %use log modulus transformation
        v3 = logmod(extraFluxes.v,10);
        v3l = logmod(extraFluxes.lb,10);
        v3l = v3 - v3l;
        v3u = logmod(extraFluxes.ub,10);
        v3u = v3u - v3;
    elseif isvar(extraFluxes,'mean') && isvar(extraFluxes,'SD')
        %use log modulus transformation
        v3 = logmod(extraFluxes.mean,10);
        v3l = logmod(extraFluxes.mean - extraFluxes.SD,10);
        v3l = v3 - v3l;
        v3u = logmod(extraFluxes.mean + extraFluxes.SD,10);
        v3u = v3u - v3;
    else
        error('Unrecognised input: extraFluxes')
    end

    %limits
    if param.measuredFluxes
        min_v = min([v;v2;v3;x-xl]);
        max_v = max([v;v2;v3;x+xu]);
    else
        min_v = min([v;v2;v3]);
        max_v = max([v;v2;v3]);
    end
    bool=~isnan(v) & ~isnan(v2) & ~isnan(v3);
end

%figure('Renderer', 'painters', 'Position', [10 10 900 1500])
%figure('units','normalized','outerposition',[0 0 1 1])
figure('units','normalized','Position',[0 0 0.5 1.5])
hold on

%patch in the background
X = [min_v max_v max_v min_v];

if param.expOrder
    Y = [nnz(negbool & bool)+nnz(zerobool & bool)+0.5 nnz(negbool & bool)+nnz(zerobool & bool)+0.5 length(measuredFluxes.mean(bool))+0.5 length(measuredFluxes.mean(bool))+0.5];
else
    Y = [nnz(negbool & bool)+nnz(zerobool & bool)+0.5 nnz(negbool & bool)+nnz(zerobool & bool)+0.5 length(predictedFluxes.v(bool))+0.5 length(predictedFluxes.v(bool))+0.5];
end

patch(X,Y,[221/255, 192/255, 237/255],'FaceAlpha',0.5)

Y = [nnz(negbool & bool)+0.5 nnz(negbool & bool)+0.5 nnz(negbool & bool)+nnz(zerobool & bool)+0.5 nnz(negbool & bool)+nnz(zerobool & bool)+0.5];
patch(X,Y,[192/255, 237/255, 207/255],'FaceAlpha',0.5)

Y = [0 0 nnz(negbool & bool)+0.5 nnz(negbool & bool)+0.5];
patch(X,Y,[192/255, 233/255, 237/255],'FaceAlpha',0.9)

if param.measuredFluxes
    if param.measuredBounds
        %errorbar
        errorbar(x(bool),(1:length(x(bool)))+0.2,xl(bool),xu(bool),'horizontal','*','LineWidth',2,'MarkerSize',12)
    else
        plot(x(bool),(1:length(x(bool)))+0.2,'*k','MarkerSize',12)%,'MarkerFaceColor',[1 1 1])
    end
end

if param.predictedBounds
    %predicted fluxes
    errorbar(v(bool),(1:length(v(bool))),vl(bool),vu(bool),'horizontal','dk','LineWidth',1.5,'MarkerSize',10)%,'MarkerFaceColor',[1 1 1])
else
    %predicted fluxes
    plot(v(bool),1:length(v(bool)),'dk','MarkerSize',12)%,'MarkerFaceColor',[1 1 1])
end

%x & y limits

if param.expOrder
    ylim([0,length(measuredFluxes.mean(bool))+0.5])
else
    ylim([0,length(predictedFluxes.v(bool))+0.5])
end

[comparisonStats,comparison] = quanQualAcc(measuredFluxes, predictedFluxes);

%other fluxes
if ~isempty(otherFluxes)
    if param.otherBounds
        errorbar(v2(bool),(1:length(v2(bool)))-0.2,v2l(bool),v2u(bool),'horizontal','o','LineWidth',1.5,'MarkerSize',10, 'MarkerEdgeColor','red','MarkerFaceColor',[1 .6 .6])%,'MarkerFaceColor',[1 1 1])
    else
        %other fluxes
        plot(v2(bool),(1:length(v2(bool)))-0.2,'o','MarkerSize',10,'MarkerEdgeColor','red','MarkerFaceColor',[1 .6 .6])
    end

    [comparisonStats,comparison] = quanQualAcc(measuredFluxes, otherFluxes);

    xlim([min_v max_v])
    %xlim([-0.5 0.5])
else
    xlim([min_v max_v])
    %xlim([-0.5 0.5])
    %title(['Comparison of \color{blue}measured \color{black}and \color{red}predicted \color{black}net exchange reaction rates. (Objective ' solution1 ')'])
end

%extra fluxes
if ~isempty(extraFluxes)
    if param.extraBounds
        errorbar(v3(bool),(1:length(v3(bool)))-0.4,v3l(bool),v3u(bool),'horizontal','x','LineWidth',1.5,'MarkerSize',10, 'MarkerEdgeColor','magenta','MarkerFaceColor',[1 .6 .6])%,'MarkerFaceColor',[1 1 1])
    else
        %extra fluxes
        plot(v3(bool),(1:length(v3(bool)))-0.4,'x','MarkerSize',10,'MarkerEdgeColor','magenta','MarkerFaceColor',[1 .6 .6])
    end

    [comparisonStats,comparison] = quanQualAcc(measuredFluxes, extraFluxes);

    xlim([min_v max_v])
    %xlim([-0.5 0.5])
else
    xlim([min_v max_v])
    %xlim([-0.5 0.5])
    %title(['Comparison of \color{blue}measured \color{black}and \color{red}predicted \color{black}net exchange reaction rates. (Objective ' solution1 ')'])
end

if isempty(param.plotTitle)
    if 0
    if param.measuredFluxes
        if isempty(otherFluxes)
            title(['Measured versus predicted (<>)' predDescription ' net exchange reaction rates.'])
        else
            if isempty(extraFluxes)
                title(['Measured versus predicted (<>)' predDescription ' and ' otherDescription '(O) net exchange reaction rates.'])
            else
                title(['Measured versus predicted (<>)' predDescription ' and ' otherDescription '(O), and ' extraDescription '(X) net exchange reaction rates.'])
            end
        end
    else
        title(['Predicted net reaction rates: ' predDescription ' vs ' otherDescription])
    end
    end
else
    title(param.plotTitle);
end

if param.measuredFluxes
    if isempty(otherFluxes)
        legend({'measured secretion','uncertain','measured uptake',expDescription,predDescription},'location','northwest');
    else
        if isempty(extraFluxes)
            legend({'measured secretion','uncertain','measured uptake',expDescription,predDescription,otherDescription},'location','northwest');
        else
            legend({'measured secretion','uncertain','measured uptake',expDescription,predDescription,otherDescription,extraDescription},'location','northwest');
        end
    end
else
    if isempty(otherFluxes)
        legend({predDescription},'location','northwest');
    else
        if isempty(extraFluxes)
            legend(['measured secretion','uncertain','measured uptake',predDescription,otherDescription],'location','northwest');
        else
            legend(['measured secretion','uncertain','measured uptake',predDescription,otherDescription,extraDescription],'location','northwest');
        end
    end
end
xlabel('Net flux (uMol/gDW/hr) on $\log_{10}$ scale: $\textrm{sign}(x) \times \log_{10}(1+x)$','interpreter','latex');

%labels
ax = gca;

%xtick=[-100:10:100];
%xticklab = cellstr(num2str(sign(xtick).*log10(1+abs(xtick)), '%d'));
%set(gca,'XTick',xtick,'XTickLabel',xticklab,'TickLabelInterpreter','tex')

if param.expOrder
    ax.YTick = 1:length(measuredFluxes.mean(bool));
    ax.YTickLabel = measuredFluxes.labels(bool);
else
    ax.YTick = 1:length(predictedFluxes.v(bool));
    ax.YTickLabel = predictedFluxes.labels(bool);
end

%hack to fix some metabolite lables
ax.YTickLabel = strrep(ax.YTickLabel,'Acetyl Isoleucine (Chemspider Id: 9964364)','Acetyl Isoleucine');
ax.YTickLabel = strrep(ax.YTickLabel,'Exchange of ','');
ax.YTickLabel = strrep(ax.YTickLabel,'3, 4-dihydroxy','3,4-dihydroxy');
ax.YTickLabel = strrep(ax.YTickLabel,'N-Acetyl-Tyrosine','N-Acetyl-L-tyrosine');
ax.YTickLabel = strrep(ax.YTickLabel,'Oxugen','Oxygen');
ax.YTickLabel = strrep(ax.YTickLabel,'Excahnge of ','');
ax.YTickLabel = strrep(ax.YTickLabel,'CYOR_u10mi','Oxidative phosphorylation (CYOR_u10mi)');
ax.YTickLabel = strrep(ax.YTickLabel,'NADH2_u10mi','NADH Dehydrogenase, Mitochondrial (NADH2_u10m)');

ax.TickLabelInterpreter='none';
ax.FontSize = 14;
switch param.labelType
    case 'metabolitePlatform'
        ax.YAxisLocation = 'right';
    case 'metabolite'
        ax.YAxisLocation = 'right';
    case 'platform'
        ax.YAxisLocation = 'left';
end

%         %right labels
%         plot(v,1:nnz(predictedExLOCB),'*r')
%         axL = gca;
%         axL.YAxisLocation = 'right';
%         axL.Color = 'none';
%         axL.Box = 'off';
%         yticks(1:length(measuredFluxes.mean))
%         axL.TickLabelInterpreter='latex';
%         axL.YColor ='black';
%         axL.FontSize = 10;
%         yt = yticks;
%         yticklabels(measuredFluxes.labelsR)

grid on
ax = gca;
ax.GridColor = [0 .5 .5];
ax.GridLineStyle = '--';
ax.GridAlpha = 0.5;
ax.Layer = 'top';

if 0
    xlim([-0.25 0.15])
    ylim([15.5 69.5])
    zoomed = 'zoom';
else
    zoomed = '';
end
if isempty(otherFluxes)
    descriptions = predDescription;
else
    descriptions = [predDescription '_vs_' otherDescription];
end

%%%% TODO %%% take the model name from the model structure
%     savefig(['iDN1_measured_vs_predicted_exchanges_' names '_' labelType zoomed])
%     saveas(gcf,['iDN1_measured_vs_predicted_exchanges_' names '_' labelType zoomed],'png')
if param.saveFigures
    savefig(['measured_vs_predicted_exchanges_' names '_' param.labelType zoomed])
    saveas(gcf,['measured_vs_predicted_exchanges_' names '_' param.labelType zoomed],'png')
end

author: opencobra

repo: cobratoolbox


text:
function [lb,ub,cand_rows] = findulb(F_struc,K)
%FINDULB Internal function to extract upper and lower variable bounds

% special code for the interval data case, to avoid overhead in the default
% function findulb

n = size(F_struc,2)-1;
lb = -inf*ones(n,1);
ub = inf*ones(n,1);
cand_rows = [];

if (K.f ~=0)
    A = -F_struc(1:K.f,2:end);
    b = F_struc(1:K.f,1);
    n = size(F_struc,2)-1;
    cand_rows = find(sum(A~=0,2)==1);
    for i = 1:length(cand_rows)
        j = find(A(cand_rows(i),:));
        ub(j)=min(ub(j),b(cand_rows(i))/A(cand_rows(i),j));
        lb(j)=max(lb(j),b(cand_rows(i))/A(cand_rows(i),j));
    end
end

if (K.l ~=0)
    A = -F_struc(K.f+1:K.f+K.l,2:end);
    b = F_struc(K.f+1:K.f+K.l,1);
    n = size(F_struc,2)-1;
    cand_rows = find(sum(A~=0,2)==1);
    for i = 1:length(cand_rows)
        j = find(A(cand_rows(i),:));
        if A(cand_rows(i),j)>0
            ub(j)=min(ub(j),sup(b(cand_rows(i))/A(cand_rows(i),j)));
        else
            lb(j)=max(lb(j),inf_(b(cand_rows(i))/A(cand_rows(i),j)));
        end
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [cmc,ind]=get_cmc_multi_cam(Y_gallery,Y_cam_gallery,Y_probe,Y_cam_probe,dist)
[~, ind]=sort(dist,2);
Y_result=Y_gallery(ind);
Y_cam_result=Y_cam_gallery(ind);
valid_probe_sample_count=0;
gallery_unique_count=length(unique(Y_gallery));
match_counter=zeros(1,gallery_unique_count);
probe_sample_count=length(Y_probe);
for i=1:probe_sample_count
    % remove gallery samples from the same camera of the probe
    Y_result_i=Y_result(i,:);
    Y_result_i(Y_cam_result(i,:)==Y_cam_probe(i))=[];
    % remove duplicated id
    Y_result_unique_i=unique(Y_result_i,'stable');
    % match for probe i
    match_i=(Y_result_unique_i==Y_probe(i));
    if sum(match_i)~=0 % if there is true matching in gallery
        valid_probe_sample_count=valid_probe_sample_count+1;
        for r=1:length(match_i)
            match_counter(r)=match_counter(r)+match_i(r);
        end
    end
end
rankk=match_counter/valid_probe_sample_count;
cmc=cumsum(rankk);
end

author: wuancong


text:
function out = convvel( in, uin, uout)
%CONVVEL Convert from velocity units to desired velocity units.
%Allowable UI and UO strings:
% 'keas' Knots equivalent airspeed
% 'm/s' metres per second
% 'kt' knots

if ~isfloat( in )
    error('Input is not floating point');
end

if nargin < 3
    uout = 'm/s';
end
    uin = lower(uin);
    uout = lower(uout);

%conversion to m/s
if strcmp('kts',uin)
    slope = 0.514444444; %m/s/kts
elseif strcmp('m/s',uin)
    slope = 1;
elseif strcmp('ft/min',uin)
    slope = 5.08e-3;
else
    error('invalid input unit string')
end

%conversion from m/s
if strcmp('kts',uout)
    slope = slope/0.514444444; %keas/m/s
elseif strcmp('ft/min',uout)
    slope = slope/5.08e-3;
%else
 %  error('invalid output unit string')
end

out = in.*slope;

author: Sable

repo: mcbench-benchmarks


text:
% Euclidean projection on a halfspace
% Sec. 8.1.1, Boyd & Vandenberghe "Convex Optimization"
% Joelle Skaf - 10/04/05
%
% The projection of x0 on a halfspace C = {x | a'*x <= b} is given by
%           minimize || x - x0 ||^2
%               s.t.    a'*x <= b
% It is also given by P_C(x0) = x0 + (b - a'*x0)*a/||a||^2 if a'*x0 > b
%                           and x0                         if a'*x0 <=b

% Input data
randn('seed',0);
n  = 10;
a  = randn(n,1);
b  = randn(1);
x0 = randn(n,1);    % a'*x0 <=b
x1 = x0 + a;        % a'*x1 > b

% Analytical solution
fprintf(1,'Computing the analytical solution for the case where a^T*x0 <=b...');
pc_x0 = x0;
fprintf(1,'Done! \n');
fprintf(1,'Computing the analytical solution for the case where a^T*x0 > b...');
pc_x1 = x1 + (b - a'*x1)*a/norm(a)^2;
fprintf(1,'Done! \n');

% Solution via QP
fprintf(1,'Computing the solution of the QP for the case where a^T*x0 <=b...');
cvx_begin quiet
    variable xs0(n)
    minimize ( norm(xs0 - x0).^2 )
    a'*xs0 <= b; %#ok
cvx_end
fprintf(1,'Done! \n');

fprintf(1,'Computing the solution of the QP for the case where a^T*x0 > b...');
cvx_begin quiet
    variable xs1(n)
    minimize ( norm(xs1 - x1).^2 )
    a'*xs1 <= b; %#ok
cvx_end
fprintf(1,'Done! \n');

% Verification
disp('-----------------------------------------------------------------');
disp('Verifying that p_C(x0) and x0_star are equal in the case where a^T*x0 <=b');
disp(['||p_C(x0) - x0_star|| = ' num2str(norm(xs0 - pc_x0))]);
disp('They should be equal to working precision');
disp('Verifying that p_C(x1) and x1_star are equal in the case where a^T*x1 > b');
disp(['||p_C(x1) - x1_star|| = ' num2str(norm(xs1 - pc_x1))]);
disp('They should be equal to working precision');

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% EUROPEAN OPTION PRICER (RUN SCRIPT)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Descritpion: Script to Price European options in Regime Switching Diffusion Models
%              using the PROJ method
% Author:      Justin Kirkby
% References:  (1) Efficient Option Pricing By Frame Duality with The Fast
%              Fourier Transform, SIAM J. Financial Math., 2015
%              (2) A unified approach to Bermudan and Barrier options under stochastic
%               volatility models with jumps. J. Economic Dynamics and Control, 2017
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[folder, name, ext] = fileparts(which( mfilename('fullpath')));
cd(folder);
addpath('../')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Step 1) CHOOSE CONTRACT/GENERAL PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
call = 1;    %For call use 1 (else, its a put)
S_0  = 100;  %Initial price
W    = 100;  %Strike            %NOTE: no error handling in place for extreme values of W (increase grid if strike falls outside)
r    = .05;  %Interest rate
q    = .00;  %dividend yield
T    = 1;    %Time (in years)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Step 2) CHOOSE MODEL PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Transition Matrix (dictates how the regimes transition)
Q = [-1 0.5 0.5;
    0.5 -1 0.5;
    0.5 0.5 -1];

sigma_vec = [0.15  0.25  0.35]; % Volatility in each state

initial_state = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Step 3) CHOOSE PROJ PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
order = 3;  %Choose spline order from { 0,1,2,3} => {Haar, Linear, Quadratic, Cubic}, Haar is least accurate
logN  = 8;   %Uses N = 2^logN  gridpoint
L1 = 8;  % determines grid witdth (usually set L1 = 8 to 15 for Levy)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% PRICE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha = L1*sqrt(T)*max(sigma_vec);   % Choose grid width based on the largest volatility
N = 2^logN;    % grid roughly centered on [- alph, alph]

tic
price = PROJ_RegimeSwitching_European(order, N, alpha, r, q, T, S_0, W, call, Q, sigma_vec, initial_state);
toc

fprintf('%.8f \n', price)


text:
report_this_filefun(mfilename('fullpath'));

%l = isnan(tmap);
%tmap(l) = 1;

%l = tmap < 150;
%tmap(l) = 150;

[lat,lon] = meshgrat(tmap,tmapleg);
%[smap,smapleg] = country2mtx('switzerland',100);
%[lat0, lon0] = meshgrat(smap,smapleg);

% tmap = km2deg(tmap/1);
[X , Y]  = meshgrid(gx,gy);

%sw = interp2(lon0,lat0,smap,lon,lat);

%  re4 = re3;

l = re4 < -6.9;
re4(l) = -6.9;
l = re4 > 6.9;
re4(l) = 6.9;

ren = interp2(X,Y,re4,lon,lat);

mi = min(min(ren));
l =  isnan(ren);
ren(l) = mi-20;

figure_w_normalized_uicontrolunits('pos',[150 500 1000 700])

hold on; axis off
axesm('MapProjection','eqaconic','MapParallels',[],...
    'MapLatLimit',[33.0 36.2],'MapLonLimit',[-118.1 -115.3])

meshm(ren,tmapleg,size(tmap),tmap);

daspectm('m',5);
tightmap
%view([10 30])
hl = camlight ; lighting phong
set(gca,'projection','perspective');

aa = a;

cd /alaskah2/stefan/ZMAP/agu99
load hec.mat

for i = 1:length(a)
    dep = interp2(lon,lat,tmap,a(i,1),a(i,2));

    pl =plot3m(a(i,2),a(i,1),dep+55,'ok');
    hold on
    fac = 64/max(a.Depth);

    facm = 6/max(a.Magnitude);
    sm = a(i,6)* facm;
    if sm < 1; sm = 1; end

    co = ceil(a(i,7)*fac)+1; if co > 63; co = 63; end
    set(pl,'Markersize',sm,'markerfacecolor','w');
end

load land.mat

for i = 1:length(a)
    dep = interp2(lon,lat,tmap,a(i,1),a(i,2));

    pl =plot3m(a(i,2),a(i,1),dep+55,'sk');
    ;
    hold on
    fac = 64/max(a.Depth);

    facm = 6/max(a.Magnitude);
    sm = a(i,6)* facm;
    if sm < 1; sm = 1; end

    co = ceil(a(i,7)*fac)+1; if co > 63; co = 63; end
    set(pl,'Markersize',sm,'markerfacecolor','w');
end
load bb.mat

for i = 1:length(a)
    dep = interp2(lon,lat,tmap,a(i,1),a(i,2));

    pl =plot3m(a(i,2),a(i,1),dep+55,'sk');

    hold on
    fac = 64/max(a.Depth);

    facm = 6/max(a.Magnitude);
    sm = a(i,6)* facm;
    if sm < 1; sm = 1; end

    co = ceil(a(i,7)*fac)+1; if co > 63; co = 63; end
    set(pl,'Markersize',sm,'markerfacecolor','w');
end

for i = 1:length(maepi(:,3))
    pl =plotm(maepi(i,2),maepi(i,1),'hk');
    set(pl,'Markersize',14,'markerfacecolor','w');
end

plotm(faults(:,2), faults(:,1),'k','Linewidth',1.0);

zdatam(handlem('allline'),10000) % keep line on surface
j = jet;
j = j(64:-1:1,:);
j = [ [ 0.85 0.9 0.9] ; j];
caxis([ min(min(re4)) max(max(re4)) ]);

colormap(j); brighten(0.1);
caxis([-7.2 7.2]);
axis off; set(gcf,'color','k')

setm(gca,'ffacecolor','k')
setm(gca,'fedgecolor','w','flinewidth',3);

setm(gca,'mlabellocation',1)
setm(gca,'meridianlabel','on')
setm(gca,'plabellocation',1)
setm(gca,'parallellabel','on')
setm(gca,'Fontcolor','w','Fontweight','bold','FontSize',12,'Labelunits','dm')

h5 = colorbar;
set(h5,'position',[0.67 0.15 0.01 0.2],'TickDir','out','Ycolor','w','Xcolor','w',...
    'Fontweight','bold','FontSize',12);
set(gcf,'Inverthardcopy','off');

a = aa;

author: CelsoReyes


text:
function [W,distances] = compute_relation(data,scale_sig,order)
%
%      [W,distances] = compute_relation(data,scale_sig)
%       Input: data= Feature_dimension x Num_data
%       ouput: W = pair-wise data similarity matrix
%              Dist = pair-wise Euclidean distance
%
%
% Jianbo Shi, 1997

% distances = zeros(size(data, 2), size(data, 2));
% for j = 1:length(data),
%   distances(j,:) = (sqrt((data(1,:)-data(1,j)).^2 +...
%                 (data(2,:)-data(2,j)).^2));
% end

distances = X2distances(data');

if (~exist('scale_sig')),
    scale_sig = 0.05*max(distances(:));
end

if (~exist('order')),
  order = 2;
end

tmp = (distances/scale_sig).^order;

W = exp(-tmp);

author: jwyang


text:
function disp(ds)
   %DISP Datasource disp overloaded operator

   % VERSION: 1.0 of waveform objects
   % AUTHOR: Celso Reyes (celso@gi.alaska.edu)
   % LASTUPDATE: 11/29/2008

   if numel(ds) > 1;
      disp(' ');
      disp(sprintf('[%s] %s object with Types:',...
         size2str(size(ds)), class(ds)));
      for j=1:numel(ds)
         alltypes(j) = {get(ds(j),'type')};
      end
      alltypes = unique(alltypes);
      for i=1:numel(alltypes)
         disp(sprintf('     %s',alltypes{i}));
      end

   elseif isunassigned(ds)
      disp('  unassigned datasource');
   else
      mytype = get(ds,'type');
      disp(['      type: ' upper(mytype)]);
      switch lower(mytype)
         case 'winston'
            disp(['    server: ' get(ds,'server')]);
            disp(['      port: ' num2str(get(ds,'port'))]);

         case 'irisdmcws'
            %do nothing

         case {'file','sac','antelope','seisan','obspy','miniseed'}
            filelist = getfilename(ds,[],[]);
            for n=1:numel(filelist)
               disp(sprintf('  location: %s',filelist{n}));
            end

         case {'user_defined'}
            filelist = getfilename(ds,[],[]);
            for n=1:numel(filelist)
               disp(sprintf('  location: %s',filelist{n}));
            end
            disp(sprintf(' Interpreting Function: %s',func2str(get(ds,'interpreter'))));
         otherwise
            warning('Datasource:unanticipatedType',...
               'unanticipated datasource type');
      end
   end;
end
function DispStr = size2str(sizeval)
   % helper function that changes the way we view the size
   %   from : [1 43 2 6] (numeric)  to  '1x43x2x6' (char)

   DispStr = sprintf('x%d', sizeval);
   DispStr = DispStr(2:end);
end

author: geoscience-community-codes

repo: GISMO


text:
function a = smoke ( n )

%*****************************************************************************80
%
%% SMOKE returns the SMOKE matrix.
%
%  Formula:
%
%    W = exp ( 2 * PI * sqrt ( -1 ) / N )
%
%    If ( J = I + 1 )
%      A(I,J) = 1
%    If ( I = N and J = 1 )
%      A(I,J) = 1
%    If ( I = J )
%      A(I,J) = W**I
%    Else
%      A(I,J) = 0
%
%  Example:
%
%    N = 5
%
%    w  1    0    0    0
%    0  w**2 1    0    0
%    0  0    w**3 1    0
%    0  0    0    w**4 1
%    1  0    0    0    w**5
%
%  Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    The matrix has an interesting spectrum.  The eigenvalues are
%    the N-th roots of unity times 2**(1/N).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Lothar Reichel, Lloyd Trefethen,
%    Eigenvalues and pseudo-eigenvalues of Toeplitz matrices,
%    Linear Algebra and Applications,
%    Volume 162-164, 1992, pages 153-185.
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, complex A(N,N), the smoke matrix.
%
  angle = 2.0 * pi / n;
  w = complex ( cos ( angle ), sin ( angle ) );

  for i = 1 : n
    for j = 1 : n
      if ( j + 1 == i )
        a(i,j) = 1.0;
      elseif ( i == n & j == 1 )
        a(i,j) = 1.0;
      elseif ( i == j )
        a(i,j) = w^i;
      else
        a(i,j) = 0.0;
      end
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% Copyright (C) 2016 Lagu
%% Copyright (C) 2017, 2019 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym zeros (@var{n})
%% @defmethodx @@sym zeros (@var{n}, @var{m})
%% Return a matrix whose elements are all 0.
%%
%% Example:
%% @example
%% @group
%% y = zeros (sym(3))
%%   @result{} y = (sym 3×3 matrix)
%%       ⎡0  0  0⎤
%%       ⎢       ⎥
%%       ⎢0  0  0⎥
%%       ⎢       ⎥
%%       ⎣0  0  0⎦
%% @end group
%% @end example
%%
%% @seealso{zeros, @@sym/ones, @@sym/eye}
%% @end defmethod

function y = zeros(varargin)

  % partial workaround for issue #13: delete when/if fixed properly
  if ((isa (varargin{nargin}, 'char')) && (strcmp (varargin{nargin}, 'sym')))
    varargin = varargin(1:(nargin-1));
  end

  if (isa (varargin{end}, 'char'))
    varargin = cell2nosyms (varargin);
    y = zeros (varargin{:});
    return
  end

  for i = 1:length(varargin)
    varargin{i} = sym(varargin{i});
  end
  if (length (varargin) == 1 && ~isscalar (varargin{1}))
    y = pycall_sympy__ ('return zeros(*_ins[0])', varargin{1});
  else
    y = pycall_sympy__ ('return zeros(*_ins)', varargin{:});
  end
end

%!test
%! y = zeros(sym(2));
%! x = [0 0; 0 0];
%! assert( isequal( y, sym(x)))

%!test
%! y = zeros(sym(2), 1);
%! x = [0; 0];
%! assert( isequal( y, sym(x)))

%!test
%! y = zeros(sym(1), 2);
%! x = [0 0];
%! assert( isequal( y, sym(x)))

%!test
%! y = zeros (sym([2 3]));
%! x = sym (zeros ([2 3]));
%! assert (isequal (y, x))

%% Check types:
%!assert( isa( zeros(sym(2), 'double'), 'double'))
%!assert( isa( zeros(3, sym(3), 'single') , 'single'))
%!assert( isa( zeros(3, sym(3)), 'sym'))
%!assert( isa( zeros(3, sym(3), 'sym'), 'sym'))

%!xtest
%! % Issue #13
%! assert( isa( zeros(3, 3, 'sym'), 'sym'))

repo: octsympy


text:
function pass = test_nonlinSys1_C2(pref)
% Test 2x2 system (sin/cos). This is nonlinearification of the test
%       test_linearSystem1
%
% Asgeir Birkisson, April 2014.

if ( nargin == 0 )
    pref = cheboppref;
end

tol = 1e-10;

% Smooth domain:
d = [-pi pi];
x = chebfun('x',d);
f = [ 0*x ; 0*x ];

%% Colloc2
pref.discretization = @chebcolloc2;

A = chebop(@(x,u,v) [u - diff(v,2) + u.^2; diff(u) + sin(v)],d);
A.lbc = @(u,v) u-1;
A.rbc = @(u,v) [v-1/2; diff(v)];

u12 = mldivide(A, f, pref);
u1 = u12{1}; u2 = u12{2};

% Want to check BCs as well.
bcFunLeft = A.lbc(u1,u2);
bcFunRight = chebfun(A.rbc(u1,u2));

pass(1) = norm( chebfun(A(x, u1, u2))) < tol;
pass(2) = norm(bcFunLeft(d(1))) < tol && norm(bcFunRight(d(end))) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% Copyright (C) 2014, 2016, 2019 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym flipud (@var{A})
%% Flip a symbolic matrix vertically.
%%
%% Example:
%% @example
%% @group
%% A = sym([1 2 pi; 4 5 6]);
%% flipud (A)
%%   @result{} (sym 2×3 matrix)
%%       ⎡4  5  6⎤
%%       ⎢       ⎥
%%       ⎣1  2  π⎦
%% @end group
%% @end example
%%
%% @seealso{@@sym/fliplr, @@sym/reshape}
%% @end defmethod

function B = flipud (A)

  cmd = { 'A, = _ins'
          'if A is None or not A.is_Matrix:'
          '    A = sp.Matrix([A])'
          'return A[::-1, :]' };

  B = pycall_sympy__ (cmd, sym(A));

end

%!test
%! % simple
%! syms x
%! A = [x 2; sym(pi) x];
%! B = [sym(pi) x; x 2];
%! assert (isequal (flipud(A), B))

%!test
%! % simple, odd # rows
%! syms x
%! A = [x 2; sym(pi) x; [1 2]];
%! B = [[1 2]; sym(pi) x; x 2];
%! assert (isequal (flipud(A), B))

%!test
%! % scalar
%! syms x
%! assert (isequal (flipud(x), x))


text:
% @author: Maziar Raissi

function params = Fractional_Levy()

clc; close all;

plt = 1;
plt_pred = 0;
save_plt = 1;

addpath ..
addpath ../Utilities
addpath ../Kernels/Fractional_Levy
addpath ../Utilities/export_fig

function CleanupFun()
    rmpath ..
    rmpath ../Utilities
    rmpath ../Kernels/Fractional_Levy
    rmpath ../Utilities/export_fig
end

finishup = onCleanup(@() CleanupFun());

rng('default')

set(0,'defaulttextinterpreter','latex')

%% Load Data
[t_star, x_star, u_star, pos] = gen_data(10^6, 0.01, 5, 100);
% u_star --> 100x5
% t_star --> 5x1
% x_star --> 100x1
N_star = size(x_star,1);
nsteps = size(t_star,1)-1;

%% Setup
N0 = 100;
N1 = 100;
%% Clean Data
i = 4;
dt = t_star(i+1) - t_star(i);

idx0 = randsample(N_star, N0);
x0 = x_star(idx0,:);
u0 = u_star(idx0,i);

idx1 = randsample(N_star,N1);
x1 = x_star(idx1,:);
u1 = u_star(idx1,i+1);

hyp = [log([1.0 0.1]) 1.1 0.0];
model = HPM(x1, u1, x0, u0, dt, hyp);
model = model.train(1500);

hyp = model.hyp;
params = hyp(3);

[pred_n_star, var_n_star] = model.predict(x_star);
var_n_star = abs(diag(var_n_star));

error = norm(pred_n_star - u_star(:,i+1))/norm(u_star(:,i+1));

fprintf(1,'=========================\n');
fprintf(1,'Step: %d, Time = %.2f\n\nNLML = %.2f, Error = %.2e\n\n', i, ...
    t_star(i+1), model.NLML, error);

str = sprintf('%.2f  ', params);
fprintf('Parameters: %s\n\n', str)
fprintf(1,'=========================\n\n');

if plt_pred == 1
    figure();
    plot_prediction_1D(x_star, u_star(:,i+1), pred_n_star, var_n_star, ...
        '$x$', '$u(t,x)$', 'Prediction');

    drawnow;
end

%% Plot Results

if plt == 1
    fig = figure();
    set(fig,'units','normalized','outerposition',[0 0 1 .5])
    subplot(3,2,1:2)
    plot(pos, 'b', 'LineWidth', 2);
    xlabel('Time')
    ylabel('Position')
    axis tight
    set(gca,'FontSize',14);
    set(gcf, 'Color', 'w');

    subplot(3,2,3);
    tit = sprintf('$t = $ %.2f\n%d training data\n', t_star(i), N0);
    bar(x_star, u_star(:,i),'m');
    xlabel('$x$')
    ylabel('$u(t,x)$')
    axis tight
    set(gca,'FontSize',14);
    set(gcf, 'Color', 'w');
    title(tit);

    subplot(3,2,4);
    tit = sprintf('$t = $ %.2f\n%d training data\n', t_star(i+1), N1);
    bar(x_star, u_star(:,i+1),'m');
    xlabel('$x$')
    ylabel('$u(t,x)$')
    axis tight
    set(gca,'FontSize',14);
    set(gcf, 'Color', 'w');
    title(tit);

    subplot(3,2,5:6);

    s = '$\begin{tabular}{|c|c|}';
    s = strcat(s, ' \hline');
    s = strcat(s, ' Correct PDE & $u_t + (-\nabla^{\sqrt{2}}_x) u = 0$ \\');
    s = strcat(s, ' \hline');
    s = strcat(s, ' Identified PDE & $u_t + (-\nabla^{', sprintf('%.3f',params(1)), '}_x) u = 0$ \\');
    s = strcat(s, ' \hline');
    s = strcat(s, ' \end{tabular}$');
    text(0.25,0.8,s,'interpreter','latex','FontSize',18)
    axis off

    if save_plt == 1
        export_fig ../Figures/Fractional_Levy.png -r300
    end

    drawnow();
end

end

function [t_star, x_star, u_star, pos] = gen_data(N,dt,m,n)

    alpha = sqrt(2);
    r = stblrnd(alpha, 0, 1, 0, N, 1);
    pos = cumsum(dt^(1/alpha)*r);

    M = 0.35;

    P = zeros(length(pos)-m,m);
    for i = 1:length(pos)-m
        y = pos(i+1:i+m) - pos(i);
        P(i,:) = y;
    end

    bins = linspace(-M,M,n+1);
    x_star = linspace(M*(1/n-1), M*(1-1/n), n)';
    t_star = dt*(1:1:m)';
    u_star = zeros(n,m);

    f = figure();
    clf
    hold
    for i = 1:m
        h = histogram(P(:,i), bins, 'Normalization', 'pdf');
        u_star(:,i) = h.Values;
    end

    close(f)

end

author: maziarraissi

repo: HPM

save_path: github-repos/MATLAB/maziarraissi-HPM


text:
function y = cvx_basis( x )
y = x.basis_;

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function [ X ] = special_add_right_not( X,S )

    s_theta=S(1:3);
    s_p=S(4:6);

    sizeS=size(S,1);
    NumberOfLandmarks=(sizeS-6)/3;

    Exps=so3_exp(s_theta);

    X.position= Exps*X.position+jaco_r(-s_theta)*s_p;

    if NumberOfLandmarks>=1
        s_landmarksMatrix=reshape(S(7:end),3,NumberOfLandmarks);
        X.landmarks(1:3,:)=Exps*X.landmarks(1:3,:)+s_landmarksMatrix;
    end

    X.orientation=Exps*X.orientation;

end

author: RomaTeng

repo: EKF-SLAM-on-Manifold


text:
classdef SimpAllInterpolationImplicit < MaterialInterpolation

    properties (Access = protected)
       dmu0
       dmu1
       dk0
       dk1
    end

    methods (Access = protected)

        function [mS,dmS] = computeMuSymbolicFunctionAndDerivative(obj)
            s.f0     = obj.matProp.mu0;
            s.f1     = obj.matProp.mu1;
            s.df0    = obj.dmu0;
            s.df1    = obj.dmu1;
            [mS,dmS] = obj.computeParameterInterpolationAndDerivative(s);
        end

        function [kS,dkS] = computeKappaSymbolicFunctionAndDerivative(obj)
            s.f0     = obj.matProp.kappa0;
            s.f1     = obj.matProp.kappa1;
            s.df0    = obj.dk0;
            s.df1    = obj.dk1;
            [kS,dkS] = obj.computeParameterInterpolationAndDerivative(s);
        end

        function [f,df] = computeParameterInterpolationAndDerivative(obj,s)
            c     = obj.computeCoefficients(s);
            rho   = sym('rho','positive');
            fSym  = obj.rationalFunction(c,rho);
            dfSym = obj.rationalFunctionDerivative(c,rho);
            f     = simplify(fSym);
            df    = simplify(dfSym);
        end

        function c = computeCoefficients(obj,s)
            f1    = s.f1;
            f0    = s.f0;
            df1   = s.df1;
            df0   = s.df0;
            c1    = sym('c1','real');
            c2    = sym('c2','real');
            c3    = sym('c3','real');
            c4    = sym('c4','real');
            coef  = [c1 c2 c3 c4];
            r1    = obj.matProp.rho1;
            r0    = obj.matProp.rho0;
            eq(1) = obj.rationalFunction(coef,r1) - f1;
            eq(2) = obj.rationalFunction(coef,r0) - f0;
            eq(3) = obj.rationalFunctionDerivative(coef,r1) - df1;
            eq(4) = obj.rationalFunctionDerivative(coef,r0) - df0;
            c     = solve(eq,[c1,c2,c3,c4]);
            c     = struct2cell(c);
            c     = [c{:}];
        end

    end

    methods (Access = protected, Static)

        function r = rationalFunction(coef,rho)
            c1  = coef(1);
            c2  = coef(2);
            c3  = coef(3);
            c4  = coef(4);
            num = (c1*rho^2 + c2*rho + 1);
            den = (c4 + rho*c3);
            r   = num/den;
        end

        function dr = rationalFunctionDerivative(coef,rho)
            c1  = coef(1);
            c2  = coef(2);
            c3  = coef(3);
            c4  = coef(4);
            n1  = c2 + 2*rho*c1;
            d1  = c4 + rho*c3;
            dr1 = n1/d1;
            n2  = -c3*(c1*rho^2 + c2*rho + 1);
            d2  = (c4 + rho*c3)^2;
            dr2 = n2/d2;
            dr  = dr1 + dr2;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


function [phi] = tdoa(arraySignals, fs, d, c)

x1 = arraySignals(:,1);
x2 = arraySignals(:,2);

[acor,lag] = xcorr(x2,x1);
[~,I] = max(abs(acor));
t = lag(I);

phi = acos(t/fs*c/d) * 180 / pi;

end

function [phi] = gcc_phat(arraySpectrum, fs, d, c)

X1 = arraySpectrum(:,1);
X2 = arraySpectrum(:,2);

G = X1 .* conj(X2);
corr = ifft(G ./ abs(G));

[~,I] = max(abs(corr));
tdoa = lag(I);

phi = acos(tdoa/fs*c/d) * 180 / pi;

end


text:
function J=polAzEquiDistRefChangeCrossGrad(pAzPt,latLonRefOld,latLonRefNew,rE)
%%POLAZEQUIDISTREFCHANGECROSSGRAD Given a point in polar azimuthal
%      coordinates with respect to a reference position on a spherical
%      Earth, find the gradient of the point in a polar azimuthal
%      coordinate system with a difference reference point taken with
%      respect to the point in the original system. This is the gradient of
%      the polAzEquidistRefChange function.
%
%INPUTS: pAzPts A 2XN set of the [ground distance; heading] points, with
%              the heading given in radians East of North, to convert.
%              Alternatively, if heights are also given, this can be a 3XN
%              set of points with the height being the third dimension.
% latLonRefOld A 2X1 [latitude;longitude] reference point in radians about
%              which pAzPts was computed.
% latLonRefNew A 2X1 [latitude;longitude] reference point in radians with
%              respect to which the function to be differentied is defined.
%           rE The radius of the reference sphere. If this argument is
%              omitted or an empty matrix is passed, the value in
%              Constants.WGS84MeanRadius is used.
%
%OUTPUTS: J A 2X2XN matrix holding the gradient of the transformation
%           evaluated at each point in pAzPts (or a 3X3XN matrix if height
%           is also provided). If one says that a point in pAzPts=[p;az;h]
%           and a converted point is [p1;az1;h1], the gradients in 3D are
%           ordered:
%           [dp1dp,   dp1dAz,  dp1dh;
%            dAz1dp, dAz1dAz, dAz1dh;
%            dh1dp,   dh1dAz,  dh1dh];
%
%The gradient is the analytic gradient of the spherical Earth conversion
%given in the function polAzEquidistRefChange. Due to a singularity, the
%gradient is not valid if the target is exactly on the opposite side of the
%Earth from the first sensor.
%
%EXAMPLE:
%Given three random points, the accuracy of this function is compared to
%numeric differentiation. The relative error is on the order of what one
%might expect due to finite precision limitations.
% rE=Constants.WGS84MeanRadius;
% latLonRef1Old=[UniformD.rand(1,[-pi/2;pi/2]);
%                UniformD.rand(1,[-pi;pi])];
% latLonRefNew=[UniformD.rand(1,[-pi/2;pi/2]);
%               UniformD.rand(1,[-pi;pi])];
% latLonPt=[UniformD.rand(1,[-pi/2;pi/2]);
%           UniformD.rand(1,[-pi;pi])];
% pAzPt=ellips2PolarAzEquidistProj(latLonPt,latLonRef1Old,rE,0);
% J=polAzEquiDistRefChangeCrossGrad(pAzPt,latLonRef1Old,latLonRefNew,rE);
% f=@(x)polAzEquidistRefChange(x,latLonRef1Old,latLonRefNew,rE,0);
% JNumDiff=numDiff(pAzPt,f,2);
% RelErr=max(max(abs((J-JNumDiff)./JNumDiff)))
%
%December 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(rE))
    rE=Constants.WGS84MeanRadius;
end

numPoints=size(pAzPt,2);
if(all(latLonRefOld==latLonRefNew))
    %For the special case of no change, force the result to be exact (avoid
    %finite precision issues).
    numDim=size(pAzPt,1);%2 or 3.
    J=repmat(eye(numDim,numDim),[1,1,numPoints]);
    return;
end

azStart=greatCircleAzimuth(latLonRefOld,latLonRefNew);
azOffset1=pi/2-azStart;
lambda0=greatCircleDistance(latLonRefOld,latLonRefNew,1);

if(size(pAzPt,1)==3)
    J=zeros(3,3,numPoints);
    J(3,3,:)=1;
else
    J=zeros(2,2,numPoints);
end

%Transform into the system with both sensors on the equator.
pAzPt(2,:)=pAzPt(2,:)+azOffset1;

p=pAzPt(1,:);
az=pAzPt(2,:);

cosPre=cos(p/rE);
sinPre=sin(p/rE);

cosLambda0=cos(lambda0);
sinLambda0=sin(lambda0);
sinAz=sin(az);
cosAz=cos(az);

denom2=1-(cosPre.*cosLambda0+sinPre.*sinAz.*sinLambda0).^2;
denom1=sqrt(denom2);

dp0dp=(cosLambda0.*sinPre-cosPre.*sinAz.*sinLambda0)./denom1;
dp0dAz=-rE.*sinPre.*cosAz.*sinLambda0./denom1;

dAz0dp=cosAz.*sinLambda0./(rE.*denom2);
dAz0dAz=sinPre.*(cosLambda0.*sinPre-cosPre.*sinAz.*sinLambda0)./denom2;

J(1,1,:)=dp0dp;
J(2,1,:)=dAz0dp;
J(1,2,:)=dp0dAz;
J(2,2,:)=dAz0dAz;

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function a = ris ( n )

%*****************************************************************************80
%
%% RIS returns the RIS matrix.
%
%  Discussion:
%
%    This is sometimes called the "dingdong" matrix, invented by F N Ris.
%
%  Formula:
%
%    A(I,J) = 1 / ( 3 + 2 * N - 2 * I - 2 * J )
%
%  Example:
%
%    N = 5
%
%    1/9  1/7  1/5  1/3   1
%    1/7  1/5  1/3   1   -1
%    1/5  1/3   1   -1  -1/3
%    1/3   1   -1  -1/3 -1/5
%     1   -1  -1/3 -1/5 -1/7
%
%  Properties:
%
%    A is a Cauchy matrix.
%
%    A is a Hankel matrix: constant along anti-diagonals.
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    The eigenvalues of A cluster around PI/2 and -PI/2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    John Nash,
%    Compact Numerical Methods for Computers: Linear Algebra and
%    Function Minimisation,
%    John Wiley, 1979, page 210.
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n
      a(i,j) = 1.0 / ( 3 + 2 * n - 2 * i - 2 * j );
    end
  end

  return
end


text:
function printDecodeResults(decodeData, candidates, candScores, alignInfo, otherInfo, params, isOutput)
  batchSize = size(candScores, 2);
  startId = decodeData.startId;

  % output translations
  [maxScores, bestIndices] = max(candScores, [], 1); % stackSize * batchSize
  for ii = 1:batchSize
    bestId = bestIndices(ii);
    translation = candidates{ii}{bestId};
    assert(isempty(find(translation>params.tgtVocabSize, 1)));

    if isOutput
      if params.forceDecoder
        printSent(params.fid, otherInfo.forceDecodeOutputs(1:decodeData.tgtLens(ii)-1, ii), params.tgtVocab, '');
      else
        printSent(params.fid, translation(1:end-1), params.tgtVocab, ''); % remove <t_eos>
      end
      if params.charTgtGen
        printSentChar(params.charFid, translation(1:end-1), params.tgtVocab, '', otherInfo.rarePositions{ii}, otherInfo.rareWords{ii}); % remove <t_eos>
      end
    end

    % log
    printSrc(params.logId, decodeData, ii, params, startId+ii-1);
    printRef(params.logId, decodeData, ii, params, startId+ii-1);
    printSent(params.logId, translation, params.tgtVocab, ['  tgt ' num2str(startId+ii-1) ': ']);
    if params.charTgtGen
      printSentChar(params.logId, translation, params.tgtVocab, ['  tgt char' num2str(startId+ii-1) ': '], otherInfo.rarePositions{ii}, otherInfo.rareWords{ii});
    end
    % align
    if params.align
      alignment = alignInfo{ii}{bestId};
      printAlign(params.logId, translation, decodeData, alignment, params, ii, startId+ii-1, 1);
      if isOutput
        printSentAlign(params.alignId, translation, decodeData, alignment, ii, params);
      end
    end
    fprintf(params.logId, '  score %g\n', maxScores(ii));

    % debug
    printSrc(2, decodeData, ii, params, startId+ii-1);
    printRef(2, decodeData, ii, params, startId+ii-1);
    printSent(2, translation, params.tgtVocab, ['  tgt ' num2str(startId+ii-1) ': ']);
    if params.forceDecoder
      printSent(2, otherInfo.forceDecodeOutputs(1:decodeData.tgtLens(ii)-1, ii), params.tgtVocab, ['  fdc ' num2str(startId+ii-1) ': ']);
    end
    if params.charTgtGen
      printSentChar(2, translation, params.tgtVocab, ['  tgt char ' num2str(startId+ii-1) ': '], otherInfo.rarePositions{ii}, otherInfo.rareWords{ii});

      % check that no words are </c_s>
      assert(isempty(find(strcmp(otherInfo.rareWords{ii}, params.tgtCharVocab{params.tgtCharEos}), 1)));
    end

    % align
    if params.align
      printAlign(2, translation, decodeData, alignment, params, ii, startId+ii-1, 1);
    end
    fprintf(2, '  score %g\n', maxScores(ii));
  end
end

function printSentAlign(fid, translation, data, alignment, ii, params)
  srcLen = data.srcLens(ii);
  if params.isReverse
    alignment = srcLen-alignment;
  end

  tgtLen = length(translation);
  for i = 1:(tgtLen-1)
    srcPos = alignment(i+1);
    assert((1 <= srcPos) && (srcPos <= srcLen), sprintf('Assertion failed: srcPos = %d, srcLen = %d', srcPos, srcLen));
    fprintf(fid, '%d-%d ', srcPos-1, i-1); % base 0
  end
  fprintf(fid, '\n');
end

function printAlign(fid, translation, data, alignment, params, ii, sentId, printWords)
  mask = data.srcMask(ii, :);
  srcLen = data.srcLens(ii);

  if params.isReverse
    alignment = srcLen-alignment;
  end

  tgtLen = length(translation);
  if printWords
    src = data.srcInput(ii,mask);
    if params.isReverse
      src = src(end:-1:1);
    end
    fprintf(fid, 'pairs %d: ', sentId);
    for i = 1:(tgtLen-1)
      srcPos = alignment(i+1);
      assert((1 <= srcPos) && (srcPos <= srcLen), sprintf('Assertion failed: srcPos = %d, srcLen = %d', srcPos, srcLen));
      fprintf(fid, '%s-%s ', params.srcVocab{src(srcPos)}, params.tgtVocab{translation(i)});
    end
    fprintf(fid, '\n');
  end

  fprintf(fid, 'align %d: ', sentId);
  for i = 1:(tgtLen-1)
    srcPos = alignment(i+1);
    assert((1 <= srcPos) && (srcPos <= srcLen));
    fprintf(fid, '%d-%d ', srcPos-1, i-1); % base 0
  end
  fprintf(fid, '\n');

end

function printSrc(fid, data, ii, params, sentId)
  mask = data.srcMask(ii,:);
  src = data.srcInput(ii,mask);
  printSent(fid, src, params.srcVocab, ['  src ' num2str(sentId) ': ']);

  if params.isReverse
    printSent(fid, src(end:-1:1), params.srcVocab, [' rsrc ' num2str(sentId) ': ']);
  end
end

function printRef(fid, data, ii, params, sentId)
  mask = data.tgtMask(ii, :);
  ref = data.tgtOutput(ii,mask);
  printSent(fid, ref, params.tgtVocab, ['  ref ' num2str(sentId) ': ']);
end

author: lmthang


text:
function hdrv = hdrvclose(hdrv)
%
%        hdrv = hdrvclose(hdrv)
%
%
%        Input:
%           -hdrv: a HDR video structure
%
%        Output:
%           -hdrv: a HDR video structure
%
%        This function closes a video stream (hdrv) for reading frames
%
%     Copyright (C) 2013-17  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

if(hdrv.streamOpen == 1)
    switch hdrv.type
        case 'TYPE_HDR_VIDEO'
            if(hdrv.permission == 'w')
                if(~isempty(hdrv.streamTMO))
                    close(hdrv.streamTMO)
                end

                if(~isempty(hdrv.streamR))
                    close(hdrv.streamR)
                end
            end
    end

    hdrv.streamOpen = 0;
end

end

author: banterle


text:
function nleqopt = mpopt2nleqopt(mpopt, model, alg)
%MPOPT2NLEQOPT   Create/modify NLEQS_MASTER options struct from MPOPT.
%
%   NLEQOPT = MPOPT2NLEQOPT(MPOPT, MODEL)
%   NLEQOPT = MPOPT2NLEQOPT(MPOPT, MODEL, ALG)
%
%   Uses a MATPOWER options struct, MPOPT, to create or modify an
%   NLEQS_MASTER options struct.
%
%   Inputs (default values in parentheses):
%       MPOPT : MATPOWER options struct
%       MODEL ('NLEQ') : (optional) one of the following model types, required
%               for selection of solver in case ALG is 'DEFAULT' (solver
%               precedence for each model type list in parentheses):
%           'NLEQ' - nonlinear equation with all continuous variables
%                   (NEWTON, FSOLVE)
%       ALG ('pf.alg') : (optional) 'pf.alg' or any valid value of
%               OPT.alg for NLEQS_MASTER. The first option indicates that
%               it should be taken from MPOPT.pf.alg.
%
%   Output:
%       NLEQOPT : an options struct for use by NLEQS_MASTER and friends
%
%   See NLEQS_MASTER, MPOPTION.

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

%% set default args
if nargin < 3
    alg = '';
    if nargin < 2
        model = '';
    end
end
if isempty(model)
    model = 'NLEQ';
else
    model = upper(model);
end

%% get ALG, MAX_IT from mpopt, if necessary
switch alg
    case {'pf.alg', ''}
        alg = upper(mpopt.pf.alg);
    otherwise
        alg = upper(alg);
end

%% create NLEQS_MASTER options struct
nleqopt = struct(   'verbose',  mpopt.verbose, ...
                    'alg',      alg, ...
                    'tol',      mpopt.pf.tol  );
switch alg
    case {'DEFAULT', 'NEWTON'}
        %% set up options
        nleqopt.newton_opt.lin_solver = mpopt.pf.nr.lin_solver;
        nleqopt.max_it = mpopt.pf.nr.max_it;
    case 'FD'
        nleqopt.max_it = mpopt.pf.fd.max_it;
    case 'GS'
        nleqopt.max_it = mpopt.pf.gs.max_it;
    case 'ZG'
        nleqopt.max_it = mpopt.pf.zg.max_it;
        nleqopt.alg = 'CORE';
    case 'FSOLVE'
        %% basic optimset options needed for fmincon
%         nleqopt.fsolve_opt.Algorithm = '';
%         nleqopt.fsolve_opt.Algorithm = 'trust-region-dogleg';
%         nleqopt.fsolve_opt.Algorithm = 'trust-region';              %% for optimoptions
%         nleqopt.fsolve_opt.Algorithm = 'trust-region-reflective';   %% for optimset
%         nleqopt.fsolve_opt.Algorithm = 'levenberg-marquardt';
        nleqopt.max_it = 0;      %% use fsolve default
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function out = Pauli(data, sicd_meta)
%PAULI Pauli Decomposition
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

if (nargin>1) && isfield(sicd_meta,'ImageFormation') && ...
        isfield(sicd_meta.ImageFormation,'TxRcvPolarizationProc')
    % Try standard linear polarizations first
    HH_ind=find(strcmpi('H:H',sicd_meta.ImageFormation.TxRcvPolarizationProc));
    HV_ind=find(strcmpi('H:V',sicd_meta.ImageFormation.TxRcvPolarizationProc));
    VH_ind=find(strcmpi('V:H',sicd_meta.ImageFormation.TxRcvPolarizationProc));
    VV_ind=find(strcmpi('V:V',sicd_meta.ImageFormation.TxRcvPolarizationProc));
    % Less common circular polarization next in not linear
    if isempty([HH_ind HV_ind VH_ind VV_ind])
        pols = cellfun(@(x) split(x,':'), sicd_meta.ImageFormation.TxRcvPolarizationProc, 'UniformOutput', false);
        HH_ind = find(cellfun(@(x) strcmpi(x{1},'LHC') && ...
            (strcmpi(x{2},'LHC') || x{2}=='H'), pols));
        HV_ind = find(cellfun(@(x) strcmpi(x{1},'LHC') && ...
            (strcmpi(x{2},'RHC') || x{2}=='V'), pols));
        VH_ind = find(cellfun(@(x) strcmpi(x{1},'RHC') && ...
            (strcmpi(x{2},'LHC') || x{2}=='V'), pols));
        VV_ind = find(cellfun(@(x) strcmpi(x{1},'RHC') && ...
            (strcmpi(x{2},'RHC') || x{2}=='H'), pols));
    end
else % Make band assumptions based on order
    switch size(data,3)
        case 2 % Co/cross
            HH_ind = 1; HV_ind = 2; VH_ind = []; VV_ind = [];
        case 3 % HH/cross/VV
            HH_ind = 1; HV_ind = 2; VH_ind = []; VV_ind = 3;
        case 4 % Full
            HH_ind = 1; HV_ind = 2; VH_ind = 3; VV_ind = 4;
    end
end

out = zeros(size(data,1),size(data,2),3);
% Co-pol,cross-pol
if xor(isscalar(HH_ind), isscalar(VV_ind)) && xor(isscalar(HV_ind), isscalar(VH_ind))
    %incoherent assignment (co is magenta, cross is green)
    out(:,:,1) = data(:,:,[HH_ind VV_ind]);
    out(:,:,2) = data(:,:,[HV_ind VH_ind]);
    out(:,:,3) = data(:,:,[HH_ind VV_ind]);
% HH/VV
elseif isscalar(HH_ind) && isscalar(VV_ind) && isempty(HV_ind) && isempty(VH_ind)
    out(:,:,1) = data(:,:,HH_ind) - data(:,:,VV_ind);
    out(:,:,2) = 0;
    out(:,:,3) = data(:,:,HH_ind) + data(:,:,VV_ind);
% HH,cross-pol,VV
elseif isscalar(HH_ind) && isscalar(VV_ind) && xor(isscalar(HV_ind), isscalar(VH_ind))
    out(:,:,1) = data(:,:,HH_ind) - data(:,:,VV_ind);
    out(:,:,2) = 2*data(:,:,[HV_ind VH_ind]);
    out(:,:,3) = data(:,:,HH_ind) + data(:,:,VV_ind);
% Full-pol (HH,HV,VH,VV)
elseif isscalar(HH_ind) && isscalar(VV_ind) && isscalar(HV_ind) && isscalar(VH_ind)
    out(:,:,1) = data(:,:,HH_ind) - data(:,:,VV_ind);
    out(:,:,2) = data(:,:,HV_ind) + data(:,:,VH_ind);
    out(:,:,3) = data(:,:,HH_ind) + data(:,:,VV_ind);
end

end

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

author: ngageoint


text:
function test6
%TEST6 test cs_reach, cs_reachr, cs_lsolve, cs_usolve
%
% Example:
%   test6
% See also: testall

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

rand ('state', 0)
maxerr = 0 ;
clf
for trial = 1:201
    n = fix (100 * rand (1)) ;
    d = 0.1 * rand (1) ;
    L = tril (sprandn (n,n,d),-1) + sprand (speye (n)) ;
    b = sprandn (n,1,d) ;

    for uplo = 0:1

        if (uplo == 1)
            % solve Ux=b instead ;
            L = L' ;
        end

        x = L\b ;
        sr = 1 + cs_reachr (L,b) ;
        sz = 1 + cs_reachr (L,b) ;

        check_if_same (sr,sz) ;

        s2 = 1 + cs_reach (L,b) ;

        if (uplo == 0)
            x3 = cs_lsolve (L,b) ;
        else
            x3 = cs_usolve (L,b) ;
        end
        x3 = x3'' ;

        spy ([L b x x3])
        drawnow

        s = sort (sr) ;
        [i j xx] = find (x) ;                                               %#ok
        [i3 j3 xx3] = find (x3) ;                                           %#ok

        if (isempty (i))
            if (~isempty (s))
                i       %#ok
                s       %#ok
                error ('!') ;
            end
        elseif (any (s ~= i))
            i       %#ok
            s       %#ok
            error ('!') ;
        end

        if (isempty (i3))
            if (~isempty (s))
                i3      %#ok
                s       %#ok
                error ('!') ;
            end
        elseif (any (s ~= sort (i3)))
            s       %#ok
            i3      %#ok
            error ('!') ;
        end

        if (any (s2 ~= sr))
            s2      %#ok
            sr      %#ok
            error ('!') ;
        end

        err = norm (x-x3,1) ;
        if (err > 1e-12)
            x       %#ok
            x3      %#ok
            uplo    %#ok
            err     %#ok
            error ('!')
        end

        maxerr = max (maxerr, err) ;

    end

    drawnow
end
fprintf ('maxerr = %g\n', maxerr) ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function det = r8po_det ( n, a_lu )

%*****************************************************************************80
%
%% R8PO_DET computes the determinant of a matrix factored by R8PO_FA.
%
%  Discussion:
%
%    The R8PO storage format is appropriate for a symmetric positive definite
%    matrix and its inverse.  (The Cholesky factor of a R8PO matrix is an
%    upper triangular matrix, so it will be in R8GE storage format.)
%
%    Only the diagonal and upper triangle of the square array are used.
%    This same storage scheme is used when the matrix is factored by
%    R8PO_FA, or inverted by R8PO_INVERSE.  For clarity, the lower triangle
%    is set to zero.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real A_LU(N,N), the LU factors from R8PO_FA.
%
%    Output, real DET, the determinant of A.
%
  det = 1.0E+00;

  for i = 1 : n
    det = det * a_lu(i,i);
  end

  return
end


text:
function [] = FinalOutput(CL, CDi, CD0, y_cp, geo, WingVolume, airfoildata, output)
% calculate score and other outputs, post to the GUI

q = 1/2*geo.density*geo.V^2;

%Determine the stall angles of attack for root and tip airfoils
alpha_r = (geo.i_r + geo.alpha)*180/pi;
chord_r = geo.c_r;
thick_r = str2num(geo.root(end-1:end));  %Determine thickness of root airfoil
Re_r = geo.V*geo.density*chord_r/geo.visc;
k = cell2mat(airfoildata{5}(geo.rootindex));
alpha_stall_r = k(1) + k(2)*log10(Re_r) + k(3)*log10(Re_r)^2;

alpha_t = (geo.i_r + geo.twist + geo.alpha)*180/pi;
chord_t = geo.c_r*geo.taper;
thick_t = str2num(geo.tip(end-1:end));  %#ok<ST2NM> %Determine thickness of tip airfoil
Re_t = geo.V*geo.density*chord_t/geo.visc;
k = cell2mat(airfoildata{5}(geo.tipindex));
alpha_stall_t = k(1) + k(2)*log10(Re_t) + k(3)*log10(Re_t)^2;

eta = (2*(1:geo.ns)-1)/(2*geo.ns);
chord = chord_r + eta*(chord_t-chord_r);
alpha = alpha_r + eta*(alpha_t-alpha_r);
alpha_stall = alpha_stall_r + eta*(alpha_stall_t-alpha_stall_r);
thick = thick_r + eta*(thick_t-thick_r);

counter = 0;    %Tracks the number of panels that are stalled
percentstalled = 0; %Tracks percentage of wing that is stalled
averagechord = 0; %Tracks average chord of stalled wing
averagealphapaststall = 0; %Tracks the average number of degrees past stall
averagethickness = 0; %Tracks the average thickness of the airfoil sections past stall
averagestallalpha = 0; %Tracks the average stall angle
for i = 1:geo.ns
    if alpha(i) > alpha_stall(i)
        counter = counter + 1;
        percentstalled = percentstalled + 1;
        averagechord = averagechord + chord(i);
        averagealphapaststall = averagealphapaststall + (alpha(i) - alpha_stall(i));
        averagethickness = averagethickness + thick(i);
        averagestallalpha = averagestallalpha + alpha_stall(i);
    end
end

if counter ~= 0 %Some part of the wing is stalled; prevents divide by zero
    percentstalled = percentstalled/geo.ns;
    averagechord = averagechord/counter;
    averagealphapaststall = averagealphapaststall/counter*pi/180;  %Convert to radians because stall model was developed in radians
    averagethickness = averagethickness/counter;
    averagestallalpha = averagestallalpha/counter;
end
%Determine lift while accounting for lift after the stall angle of attack
c = averagethickness/12;   %Thickness coefficient
d = 100*c;  %Determined by trial and error to closely model the shape and behavior of the lift curve
k = 400*c;
CL_stall = sum(CL.c+CL.alpha*averagestallalpha*pi/180);  %Average max lift coefficient
%The following equation is the result of modeling the lift
%coefficient of the airfoil after the stall angle of attack as a
%2nd order ODE with initial position of CL_stall and initial
%velocity of CL_alpha.  See Stall Model.nb for Mathematica
%derivation
if percentstalled == 0
    stallcorrection = 1;
else
    stallcorrection = exp(0.5*(-d-sqrt(d^2-4*k))*(averagealphapaststall))*(-2*CL.alpha-CL_stall*d+CL_stall*sqrt(d^2-4*k))/(2*sqrt(d^2-4*k))...
        +exp(0.5*(-d+sqrt(d^2-4*k))*(averagealphapaststall))*(2*CL.alpha+CL_stall*d+CL_stall*sqrt(d^2-4*k))/(2*sqrt(d^2-4*k));
end

CL = ((1-percentstalled)*(CL.c+CL.alpha*geo.alpha)+percentstalled*stallcorrection);

L = CL*q*geo.S;

Dind = CDi*q*geo.S;

%CD0 is profile drag of wing, S_Sref*Cfe is the zero-lift drag of the
%fuselage from section 2.9.1 of Anderson, Aircraft Performance and Design
%The sum of these two terms represents the total skin friction drag on the
%wing body or the zero-lift drag
Dprofile = (CD0 + geo.S_Sref*geo.cf)*q*geo.S;

Dtotal = Dind + Dprofile;  %Newton

BendMoment = -(y_cp.c+y_cp.alpha*geo.alpha)*L;

%Determine useful payload
Fuel_weight = geo.wingfuel*geo.fueldens*WingVolume*9.81;   %Fuel weight in N
Payload = L - geo.emptyweight - Fuel_weight;    %Useful payload in N

%Range equations from Aircraft Performance and Design
switch geo.engine
    case 'prop'
        Range = geo.propefficiency/geo.propSFC*L/Dtotal*log((geo.emptyweight+Fuel_weight)/geo.emptyweight);  %25% usable volume
        SFC = geo.propSFC;
    case 'jet'
        CD = CD0 + CDi;
        Range = 2/geo.jetTSFC*sqrt(2/(geo.S*geo.density)*CL/CD^2)*(sqrt(geo.emptyweight+Fuel_weight)-sqrt(geo.emptyweight));
        SFC = geo.jetTSFC;
    otherwise
        Range = 0;
        SFC = 0;
end

%Score function
Baseline_payload = 45000; %C130 maximum payload, Wikipedia
Cost_payload = 10000/400; %An extra 10000 lbf of lift is worth 1000 miles
Baseline_moment = 422000; %C130 moment from Wing Designer using C130 geometry
Cost_moment = 400000/10000; %Doubling the baseline moment incurs a penalty of 10000 miles
Baseline_airspeed = 292; %292 knots (C130 cruise), Wikipedia
Cost_airspeed = 10/100; %An extra 10 knots of cruising speed is worth 100 miles
Baseline_span = 132/3.28; %C130 wing span, Wikipedia
Cost_span = 20/3.28/2000; %An extra 20 feet incurs a penalty of 1200 miles
Score = Range/1609 + (Payload*0.22481-Baseline_payload)/Cost_payload - (BendMoment - Baseline_moment)/Cost_moment + ...
    (geo.V/0.5144 - Baseline_airspeed)/Cost_airspeed - (geo.b - Baseline_span)/Cost_span; %Convert range to miles; lift to lbf; airspeed to knots

	%Set output
set(output.bendmoment,'String',num2str(round(BendMoment)));
set(output.totaldrag,'String',num2str(round(Dtotal*0.22481*10)/10));%Convert to lbf
set(output.inddrag,'String',num2str(round(Dind*0.22481*10)/10));%Convert to lbf
set(output.profdrag,'String',num2str(round(Dprofile*0.22481*10)/10));%Convert to lbf
set(output.range,'String',num2str(round(Range/1609*10)/10));%Convert m to miles
set(output.volume,'String',num2str(round(WingVolume*35.315*100)/100));%Convert to ft^3
set(output.calctime,'String',num2str(etime(clock,geo.t0)));
set(output.lift,'String',num2str(round(L*0.22481)));  %Convert to lbf
set(output.fuelweight,'String',num2str(round(Fuel_weight*0.22481)));  %Convert to lbf
set(output.payload,'String',num2str(round(Payload*0.22481)));  %Convert to lbf

%Conditional formatting; if payload < requirements, font -> red, score -> 0
if (Payload*0.22481 < 45000) ||  (Range/1609 < 2360) %Baseline payload and range of C130, Wikipedia
    set(output.SCORE,'String','0');
else
    set(output.SCORE,'String',num2str(round(Score*10)/10));
end

if Payload*0.22481 < 45000
    set(output.payload,'ForegroundColor','red')
else
    set(output.payload,'ForegroundColor','black')
end

if Range/1609 < 2360
    set(output.range,'ForegroundColor','red')
else
    set(output.range,'ForegroundColor','black')
end

author: Sable


text:
function varargout = sound_field_mono_wfs(X,Y,Z,xs,src,f,conf)
%SOUND_FIELD_MONO_WFS sound field for WFS
%
%   Usage: [P,x,y,z,x0] = sound_field_mono_wfs(X,Y,Z,xs,src,f,conf)
%
%   Input parameters:
%       X           - x-axis / m; single value or [xmin,xmax] or nD-array
%       Y           - y-axis / m; single value or [ymin,ymax] or nD-array
%       Z           - z-axis / m; single value or [zmin,zmax] or nD-array
%       xs          - position of virtual source / m
%       src         - source type of the virtual source
%                         'pw' - plane wave (xs is the direction of the
%                                plane wave in this case)
%                         'ps' - point source
%                         'fs' - focused source
%       f           - monochromatic frequency / Hz
%       conf        - configuration struct (see SFS_config)
%
%   Output parameters:
%       P           - simulated sound field
%       x           - corresponding x values / m
%       y           - corresponding y values / m
%       z           - corresponding z values / m
%       x0          - active secondary sources / m
%
%   SOUND_FIELD_MONO_WFS(X,Y,Z,xs,src,f,conf) simulates a monochromatic sound
%   field for the given source type (src) synthesized with wave field synthesis
%   for the frequency f.
%
%   To plot the result use:
%   plot_sound_field(P,X,Y,Z,x0,conf);
%   or simple call the function without output argument:
%   sound_field_mono_wfs(X,Y,Z,xs,src,f,conf)
%
%   See also: plot_sound_field, sound_field_imp_wfs, driving_function_mono_wfs

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 7;
nargmax = 7;
narginchk(nargmin,nargmax);
isargxs(xs);
isargpositivescalar(f);
isargchar(src);
isargstruct(conf);

%% ===== Configuration ==================================================
if strcmp('2D',conf.dimension)
    greens_function = 'ls';
else
    greens_function = 'ps';
end

%% ===== Computation ====================================================
% Get the position of the loudspeakers and its activity
x0 = secondary_source_positions(conf);
x0 = secondary_source_selection(x0,xs,src);
x0 = secondary_source_tapering(x0,conf);
% Driving function
D = driving_function_mono_wfs(x0,xs,src,f,conf);
% Wave field
[varargout{1:min(nargout,4)}] = ...
    sound_field_mono(X,Y,Z,x0,greens_function,D,f,conf);
% Return secondary sources if desired
if nargout==5, varargout{5}=x0; end

author: sfstoolbox

repo: sfs-matlab


text:
% TraceHeaderDef : Defines names, position, and precision for Trace Headers
%
% % To get a Matlab structure with trace header definitions call:
% STH==TraceHeaderDef;
% % To get a list fo trace header definision listed on the screen call:
% STH==TraceHeaderDef(1)
%
% See also: ReadSegyTraceHeaderValue, WriteSegyTraceHeaderValue
%
function STH=TraceHeaderDef(print);

if nargin>0
    STH=TraceHeaderDef;
    fn=fieldnames(STH);
    disp(sprintf('%4s %6s %s','POS','PREC','Traece Header Name'))
    for i=1:length(fn)
        disp(sprintf('%4d %6s %s',STH.(fn{i}).pos,STH.(fn{i}).precision,fn{i}))
    end
    return
end

%NEW
STH.TraceSequenceLine.pos=0; STH.TraceSequenceLine.precision='int32';
STH.TraceSequenceFile.pos=4; STH.TraceSequenceFile.precision='int32';
STH.FieldRecord.pos=8;       STH.FieldRecord.precision='int32';
STH.TraceNumber.pos=12;      STH.TraceNumber.precision='int32';
STH.EnergySourcePoint.pos=16;STH.EnergySourcePoint.precision='int32';
STH.cdp.pos=20;              STH.cdp.precision='int32';
STH.cdpTrace.pos=24;         STH.cdpTrace.precision='int32';

STH.TraceIdentifactionCode.pos=28;  STH.TraceIdentifactionCode.precision='int16';
STH.NSummedTraces.pos=30;           STH.NSummedTraces.precision='int16';
STH.NStackedTraces.pos=32;          STH.NStackedTraces.precision='int16';
STH.DataUse.pos=34;                 STH.DataUse.precision='int16';

STH.offset.pos=36;STH.offset.precision='int32';
STH.ReceiverGroupElevation.pos=40;STH.ReceiverGroupElevation.precision='int32';
STH.SourceSurfaceElevation.pos=44;STH.SourceSurfaceElevation.precision='int32';
STH.SourceDepth.pos=48;STH.SourceDepth.precision='int32';
STH.ReceiverDatumElevation.pos=52;STH.ReceiverDatumElevation.precision='int32';
STH.SourceDatumElevation.pos=56;STH.SourceDatumElevation.precision='int32';
STH.SourceWaterDepth.pos=60;STH.SourceWaterDepth.precision='int32';
STH.GroupWaterDepth.pos=64;STH.GroupWaterDepth.precision='int32';

STH.ElevationScalar.pos=68;STH.ElevationScalar.precision='int16';
STH.SourceGroupScalar.pos=70;STH.SourceGroupScalar.precision='int16';

STH.SourceX.pos=72;STH.SourceX.precision='int32';
STH.SourceY.pos=76;STH.SourceY.precision='int32';
STH.GroupX.pos=80;STH.GroupX.precision='int32';
STH.GroupY.pos=84;STH.GroupY.precision='int32';

STH.CoordinateUnits.pos=88;STH.CoordinateUnits.precision='int16';
STH.WeatheringVelocity.pos=90;STH.WeatheringVelocity.precision='int16';
STH.SubWeatheringVelocity.pos=92;STH.SubWeatheringVelocity.precision='int16';
STH.SourceUpholeTime.pos=94;STH.SourceUpholeTime.precision='int16';
STH.GroupUpholeTime.pos=96;STH.GroupUpholeTime.precision='int16';
STH.SourceStaticCorrection.pos=98;STH.SourceStaticCorrection.precision='int16';
STH.GroupStaticCorrection.pos=100;STH.GroupStaticCorrection.precision='int16';
STH.TotalStaticApplied.pos=102;STH.TotalStaticApplied.precision='int16';
STH.LagTimeA.pos=104;STH.LagTimeA.precision='int16';
STH.LagTimeB.pos=106;STH.LagTimeB.precision='int16';
STH.DelayRecordingTime.pos=108;STH.DelayRecordingTime.precision='int16';
STH.MuteTimeStart.pos=110;STH.MuteTimeStart.precision='int16';
STH.MuteTimeEND.pos=112;STH.MuteTimeEND.precision='int16';

STH.ns.pos=114;STH.ns.precision='uint16';
STH.dt.pos=116;STH.dt.precision='uint16';

STH.GainType.pos=118;STH.GainType.precision='int16';
STH.InstrumentGainConstant.pos=120;STH.InstrumentGainConstant.precision='int16';
STH.InstrumentInitialGain.pos=122;STH.InstrumentInitialGain.precision='int16';
STH.Correlated.pos=124;STH.Correlated.precision='int16';
STH.SweepFrequenceStart.pos=126;STH.SweepFrequenceStart.precision='int16';
STH.SweepFrequenceEnd.pos=128;STH.SweepFrequenceEnd.precision='int16';
STH.SweepLength.pos=130;STH.SweepLength.precision='int16';
STH.SweepType.pos=132;STH.SweepType.precision='int16';
STH.SweepTraceTaperLengthStart.pos=134;STH.SweepTraceTaperLengthStart.precision='int16';
STH.SweepTraceTaperLengthEnd.pos=136;STH.SweepTraceTaperLengthEnd.precision='int16';
STH.TaperType.pos=138;STH.TaperType.precision='int16';
STH.AliasFilterFrequency.pos=140;STH.AliasFilterFrequency.precision='int16';
STH.AliasFilterSlope.pos=142;STH.AliasFilterSlope.precision='int16';
STH.NotchFilterFrequency.pos=144;STH.NotchFilterFrequency.precision='int16';
STH.NotchFilterSlope.pos=146;STH.NotchFilterSlope.precision='int16';
STH.LowCutFrequency.pos=148;STH.LowCutFrequency.precision='int16';
STH.HighCutFrequency.pos=150;STH.HighCutFrequency.precision='int16';
STH.LowCutSlope.pos=152;STH.LowCutSlope.precision='int16';
STH.HighCutSlope.pos=154;STH.HighCutSlope.precision='int16';
STH.YearDataRecorded.pos=156;STH.YearDataRecorded.precision='int16';
STH.DayOfYear.pos=158;STH.DayOfYear.precision='int16';
STH.HourOfDay.pos=160;STH.HourOfDay.precision='int16';
STH.MinuteOfHour.pos=162;STH.MinuteOfHour.precision='int16';
STH.SecondOfMinute.pos=164;STH.SecondOfMinute.precision='int16';
STH.TimeBaseCode.pos=166;STH.TimeBaseCode.precision='int16';
%STH.TimeBaseCodeText.pos=168;STH.TimeBaseCodeText.precision='int16';
STH.TraceWeightningFactor.pos=168;STH.TraceWeightningFactor.precision='int16';
STH.GeophoneGroupNumberRoll1.pos=170;STH.GeophoneGroupNumberRoll1.precision='int16';
STH.GeophoneGroupNumberFirstTraceOrigField.pos=172;STH.GeophoneGroupNumberFirstTraceOrigField.precision='int16';
STH.GeophoneGroupNumberLastTraceOrigField.pos=174;STH.GeophoneGroupNumberLastTraceOrigField.precision='int16';
STH.GapSize.pos=176;STH.GapSize.precision='int16';
STH.OverTravel.pos=178;STH.OverTravel.precision='int16';

STH.cdpX.pos=180;STH.cdpX.precision='int32';
STH.cdpY.pos=184;STH.cdpY.precision='int32';
STH.Inline3D.pos=188;STH.Inline3D.precision='int32';
STH.Crossline3D.pos=192;STH.Crossline3D.precision='int32';
STH.ShotPoint.pos=196;STH.ShotPoint.precision='int32';

STH.ShotPointScalar.pos=200;STH.ShotPointScalar.precision='int16';

STH.TraceValueMeasurementUnit.pos=202;STH.TraceValueMeasurementUnit.precision='int16';

STH.TransductionConstantMantissa.pos=204;STH.TransductionConstantMantissa.precision='int32';

STH.TransductionConstantPower.pos=208;STH.TransductionConstantPower.precision='int16';
STH.TransductionUnit.pos=210;STH.TransductionUnit.precision='int16';
STH.TraceIdentifier.pos=212;STH.TraceIdentifier.precision='int16';
STH.ScalarTraceHeader.pos=214;STH.ScalarTraceHeader.precision='int16';
STH.SourceType.pos=216;STH.SourceType.precision='int16';

STH.SourceEnergyDirectionMantissa.pos=218;STH.SourceEnergyDirectionMantissa.precision='int32';

STH.SourceEnergyDirectionExponent.pos=222;STH.SourceEnergyDirectionExponent.precision='int16';

STH.SourceMeasurementMantissa.pos=224;STH.SourceMeasurementMantissa.precision='in32';

STH.SourceMeasurementExponent.pos=228;STH.SourceMeasurementExponent.precision='int16';
STH.SourceMeasurementUnit.pos=230;STH.SourceMeasurementUnit.precision='int16';

STH.UnassignedInt1.pos=232;STH.UnassignedInt1.precision='int32';
STH.UnassignedInt2.pos=236;STH.UnassignedInt2.precision='int32';

author: cultpenguin

repo: segymat

save_path: github-repos/MATLAB/cultpenguin-segymat


text:
function varargout=lomb(varargin)
%
% [Pxx,F]=lomb(x,dcOffset,smooth)
%
%    Wrapper to WFDB LOMB:
%         http://www.physionet.org/physiotools/wag/lomb-1.htm
%
% Transforms a real-valued time series 'x' into a power spectrum 'X', using a
% technique known as the Lomb periodogram. The input is a Nx2 matrix containing
% a sampled time series, presented as two columns of numbers (the sample times
% in the first column and the sample values in the second). The intervals between
% consecutive samples need not be uniform.
%
%Input Parameters:
% x
%       Nx2 vector of doubles. First column is sample time index (in
%       seconds), and second column is the sample value of the signal at
%       that time.
%
% dcOffset (Optional)
%       Booelan. If present add constant to input samples ( x(:,2) ), such that the mean
%       values of the time series is zero (default=1).
%
% smooth   (Optional)
%       Boolean String specifying the if the output should be smoothed (default =1).
%
%
%Output Parameters:
%
%Pxx
%       Mx1 Double. Estimated power spectrum.
%
%F
%       Mx1 Double. Frequency of the estimated power spectrum (Hz).
%
%
% CITING CREDIT: To credit this function, please cite the following paper at your work:
%
%Moody, G.B.
%    Spectral analysis of heart rate without resampling. Computers in Cardiology 1993, pp. 715-718 (IEEE Computer Society Press, 1993). http://www.physionet.org/physiotools/lomb/lomb.html .
%
%
%Additional References:
%Lomb, N.R.
%    Least-squares frequency analysis of unequally spaced data. Astrophysics and Space Science 39:447-462 (1976).
%Press, W.H, and Rybicki, G.B.
%    Fast algorithm for spectral analysis of unevenly sampled data. Astrophysical J. 338:277-280 (1989).
%Press, W.H. Teukolsky, S.A., Vetterling, W.T., and Flannery, B.P.
%    Numerical Recipes in C: the Art of Scientific Computing, pp. 575-584 (Cambridge Univ. Press, 1992).
%Moody, G.B.
%    Spectral analysis of heart rate without resampling. Computers in Cardiology 1993, pp. 715-718 (IEEE Computer Society Press, 1993). http://www.physionet.org/physiotools/lomb/lomb.html .
%
%
%
% MATLAB wrapper written by Ikaro Silva, 2013
% Last Modified: -
% Version 1.0
% Since 0.9.0
%
%
% %Example: Heart Rate Spectral Analysis:
% [tm, signal]=rdsamp('mitdb/100',1);
% [ann]=rdann('mitdb/100','atr');
% [Pxx,F]=lomb([tm(ann) signal(ann)]);
% plot(F,Pxx);grid on;hold on
%
% See also RDANN, TACH, SQRS, WQRS

%endOfHelp

persistent javaWfdbExec
if(isempty(javaWfdbExec))
    javaWfdbExec=getWfdbClass('lomb');
end

%Set default pararamter values
%[Pxx,F]
inputs={'x','dcOffset','smooth'};
dcOffset=1;
smooth=1;
for n=1:nargin
    if(~isempty(varargin{n}))
        eval([inputs{n} '=varargin{n};'])
    end
end

wfdb_argument={'-P'};

if(dcOffset)
     wfdb_argument{end+1}='-z';
end
if(smooth)
     wfdb_argument{end+1}='-s';
end

wfdb_argument{end+1}='-';
del=repmat([' '],size(x(:,1)));
data=[num2str(x(:,1)) del num2str(x(:,2))];
javaWfdbExec.setArguments(wfdb_argument);
pxx=char(javaWfdbExec.execWithStandardInput(cellstr(data)));
pxx=sscanf(pxx(2:end-1), '%f %f,');

varargout{1}=pxx(2:2:end);
varargout{2}=pxx(1:2:end);

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function x = project_box(v, l, u)
% PROJECT_BOX    Project a point onto a box (hyper-rectangle).
%
%   project_box(v,l,u) is the projection of v onto
%   the set { x | l <= x <= u }.

    x = max(l, min(v, u));
end

author: cvxgrp

repo: proximal

save_path: github-repos/MATLAB/cvxgrp-proximal


text:
classdef layerTotalAvgPool

    properties
        type= 'custom'
        name= 'avgTotal'
        precious= false
    end

    methods

        function l= layerTotalAvgPool(name)
            if nargin>0, l.name= name; end
        end

    end

    methods (Static)

        function res1= forward(p, res0, res1)
            res1.x= vl_nnpool(res0.x, [size(res0.x,1), size(res0.x,2)], ...
                'method', 'avg');
        end

        function res0= backward(p, res0, res1)
            res0.dzdx= vl_nnpool(res0.x, [size(res0.x,1), size(res0.x,2)], res1.dzdx, ...
                'method', 'avg');
        end

    end

end

author: Relja

repo: netvlad

save_path: github-repos/MATLAB/Relja-netvlad


text:
function c = demean(c)

   % DEMEAN removes the mean of each trace.
   %
   % C = DEMEAN(C) removes the DC offset from each trace. In most cases it is
   % unnecessary to call this function directly. By default, all traces are
   % demeaned and detrended when they are loaded into a correlation object.
   % This is one of the assumptions of the correlation toolbox.

   % Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   c.traces = demean(c.traces);
end

author: geoscience-community-codes


text:
function [O ,O2, DzDw1 , DzDw2 ] = xmid(I1, I2, Y, gamma, Rho, DzDy)
%% Reconstruction layer in the middle of the network.
%% This layer has two parameters: H_{l} = sum_{m=1}^{s} gamma_{l,m}B_{m};  \rho_{l}
%% This layer has two inputs: I1 = z ; I2 = \beta
%% Copyright (c) 2017 Yan Yang
%% All rights reserved.

%% network setting
config;
fN = nnconfig.FilterNumber;
fS = nnconfig.FilterSize;
pad = nnconfig.Padding;
gp = nnconfig.EnableGPU;
s = fS*fS-1;

B=filter_base( );
for i=1:fN
    H(:,:,i) = reshape( B * gamma( : ,  i) , fS ,fS);
end
HT = rot90(H,2);
load('./mask.mat')
mask = logical( ifftshift(mask) );
Denom1 = zeros(256 , 256) ; Denom1(mask) = 1 ;
[m ,n] = size(Y);
Denom2=zeros(m,m);

for k=1:fN
    prd = sqrt(Rho(k));
    Denom2 = Denom2 + abs( psf2otf (prd * H(:,:,k),[m,n])).^2;
end
Denom = Denom1+Denom2;
Denom(find(Denom == 0)) = 1e-6;
Q2=1./Denom;

if nargin == 5
    Pr=zeros(m,n);
    for i = 1:fN
        tp1 = I1(:,:,i) - I2(:,:,i);
        tp2 = imfilter(double(tp1),double(HT(:,:,i)),'same','circular','conv');
        Pr= Pr+Rho(i)*tp2;
    end
    O = real( ifft2(( fft2 ( Pr ) + Y ) .* Q2));
end

if nargin == 6

    if gp
        Hp = gpuArray(H);
        [Bj1p,PSp] = FBASE(B);
        Q2=gpuArray(Q2);
        DzDy=gpuArray(DzDy);
        H1p=gpuArray(HT);
        Y=gpuArray(Y);
        Ds = gpuArray(zeros(m,m));
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  O1 O2
        A = (-1)*( Q2 .* Q2);
        Trans1 =  ifft2(Q2.*fft2(DzDy));
        Trans1 = real(Trans1);
        Trans=gather(Trans1);
        PadT = padImage_g(Trans,[pad,pad],'circular');
        PadT = gpuArray(PadT);
        for i = 1:fN
            S = conv2(PadT,Hp(:,:,i),'valid');
            S = gather(S);
            O(:,:,i) = Rho(i)*S;
            O2(:,:,i) = (-1)*Rho(i)*S;
            tp5 = I1(:,:,i)-I2(:,:,i);
            aaa = padImage_g( tp5,[pad,pad],'circular');
            aaa = gpuArray(aaa);
            tp52 = conv2(aaa,H1p(:,:,i),'valid');
            Ds= Ds +Rho(i) *tp52;
        end
        Njj = fft2(Ds)+Y;
        for i = 1:fN
            PS1 = psf2otf(Hp(:,:,i), [m,m]);
            tp5 = I1(:,:,i)-I2(:,:,i);
            Padtp = padImage_g(tp5,[pad,pad],'circular');
            Padtp = gpuArray(Padtp);
            tp54 = conv2(Padtp,H1p(:,:,i),'valid');
            tp54 = fft2(tp54);
            Nii2 = abs(PS1).^2 ;
            tp1 = A.* Nii2 .*Njj;
            tp2 = Q2.*tp54;
            tp1 = ifft2(tp1);
            tp2 = ifft2(tp2);
            temp = DzDy.*(tp1+tp2);
            temp = real(gather(temp));
            DzDw2(i) = sum(temp(:));
            for  j = 1 : s
                PS3 = 2*(PSp(:,:,j).*PS1);
                tp6 = A.*PS3.*Njj;
                PS3 = ifft2(tp6);
                DS = conv2(Padtp,Bj1p(:,:,j),'valid');
                DS = fft2(DS);
                DS = Q2.*DS;
                DS = ifft2(DS);
                PfN = DzDy.*(PS3+DS);
                SSS = real(PfN*Rho(i));
                DzDw1(j,i) = sum(SSS(:));
            end

        end
        DzDw1=gather(DzDw1);
    else
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  O1, O2
        for k = 1:fN
            Trans =  ifft2(Q2.*fft2(DzDy));
            S = imfilter(double(Trans),double(H(:,:,k)),'same','circular','conv');
            O(:,:,k) = Rho(k)*S;
            O2(:,:,k) = (-1)*Rho(k)*S;
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DzDw1
        A=(-1)* Q2 .* Q2;
        Ds=zeros(m,m);
        for i=1:fN
            tp1 = I1(:,:,i)-I2(:,:,i);
            tp3=Rho(i) * imfilter(double(tp1) ,double(HT(:,:,i)),'same','circular','conv');
            Ds= Ds +tp3;
        end
        Njj=fft2(Ds)+Y;
        for i=1:fN
            for j=1:s
                Bj=reshape(B(:,j),fS,fS);
                Bj1=rot90(Bj,2);
                PS=psf2otf(Bj1, [m,m]);
                PS1=psf2otf(H(:,:,i), [m,m]);
                PS3=2*PS.*PS1;
                PS3=DzDy.*ifft2(Rho(i)*A.*PS3.*Njj);
                tp1=I1(:,:,i)-I2(:,:,i);
                DS= imfilter(double(tp1) ,double(Bj1),'same','circular','conv');
                DS=DzDy.*ifft2(Rho(i)*Q2.*fft2(DS));
                DzDw1(j,i)=sum(sum(PS3+DS));
            end
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DzDw2
        for k=1:fN
            Nii =abs(psf2otf(H(:,:,k),[m,n])).^2 ;
            tp1=I1(:,:,k) - I2(:,:,k);
            tp4 = imfilter(double(tp1),double(HT(:,:,k)),'same','circular','conv');
            tp4 =fft2(tp4);
            temp1=DzDy.*ifft2(A.*Nii.*Njj );
            temp2=DzDy.*ifft2( Q2.*tp4);
            temp=sum(temp1(:))+sum(temp2(:));
            DzDw2(k)=temp;
        end
    end
end
end

repo: Deep-ADMM-Net


text:
function linpack_c_test23 ( )

%*****************************************************************************80
%
%% TEST23 tests CPPCO.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2009
%
%  Author:
%
%    John Burkardt
%
  n = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST23\n' );
  fprintf ( 1, '  For a single precision complex (C)\n' );
  fprintf ( 1, '  Hermitian positive definite packed matrix (PP),\n' );
  fprintf ( 1, '  CPPCO estimates the reciprocal condition number.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The matrix order is N = %d\n', n );
%
%  Set the values of the matrix A.
%
  a(1) = complex ( 2.5281,  0.0000 );

  a(2) = complex ( 2.1341, -0.2147 );
  a(3) = complex ( 3.0371,  0.0000 );

  a(4) = complex ( 2.4187,  0.2932 );
  a(5) = complex ( 2.0905,  1.1505 );
  a(6) = complex ( 2.7638,  0.0000 );
%
%  Estimate the condition.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Estimate the condition number.\n' );

  [ a, rcond, info ] = cppco ( a, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Reciprocal condition number = %f\n', rcond );

  return
end

author: johannesgerer


text:
%%***********************************************************
%% lmiexamp1: generate SDP data for the following LMI problem
%%
%%  max  -eta
%%  s.t. B*P + P*B'  <= 0
%%        -P         <= -I
%%         P - eta*I <= 0
%%         P(1,1)     = 1
%%***********************************************************
%% Here is an example on how to use this function to
%% find an optimal P.
%%
%% B = [-1  0  0; 5 -2  0; 1  1 -1];
%% [blk,At,C,b] = lmiexamp1(B);
%% [obj,X,y,Z] = sqlp(blk,At,C,b);
%% P = smat(blk(1,:),y);
%%
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 20 Apr 02
%%***********************************************************

   function [blk,At,C,b] = lmiexamp1(B);

   n = length(B); n2 = n*(n+1)/2;
   I  = speye(n);
   z0 = sparse(n2,1);
   blktmp{1,1} = 's'; blktmp{1,2} = n;
%%
   blk{1,1} = 's'; blk{1,2} = n;
   blk{2,1} = 's'; blk{2,2} = n;
   blk{3,1} = 's'; blk{3,2} = n;
   blk{4,1} = 'u'; blk{4,2} = 1;
%%
   At{1,1} = [lmifun(B,I),     z0];
   At{2,1} = [lmifun(-I/2,I),  z0];
   At{3,1} = [lmifun(I/2,I),   svec(blktmp,-I,1)];
   At{4,1} = sparse([1, zeros(1,n2)]);
%%
   C{1,1} = sparse(n,n);
   C{2,1} = -speye(n);
   C{3,1} = sparse(n,n);
   C{4,1} = 1;
%%
   b = [zeros(n2,1); -1];
%%**********************************************************

author: zarathustr

repo: LibQPEP