text:
% DESCRIPTION: Demo code for running seqNMF on simulated and real data,
% including how to test significance of each factor on held-out data, and
% how to select lambda
%
% ------------------------------------------------------------------------
% Andrew Bahle and Emily Mackevicius 1.26.2018
%
% See paper:
% https://www.biorxiv.org/content/early/2018/03/02/273128
%% Generate some synthetic data
number_of_seqences = 3;
T = 3000; % length of data to generate
Nneurons = 10*ones(number_of_seqences,1); % number of neurons in each sequence
Dt = 3.*ones(number_of_seqences,1); % gap between each member of the sequence
NeuronNoise = 0.001; % probability of added noise in each bin
SeqNoiseTime = zeros(number_of_seqences,1); % Jitter parameter = 0%
SeqNoiseNeuron = 1.*ones(number_of_seqences,1); % Participation parameter = 100%
X = generate_data(T,Nneurons,Dt,NeuronNoise,SeqNoiseTime,SeqNoiseNeuron,0,0,0,0,0);

%% Fit with seqNMF
K = 5;
L = 50;
lambda =.005;
shg; clf
display('Running seqNMF on simulated data (2 simulated sequences + noise)')
[W,H] = seqNMF(X,'K',K, 'L', L,'lambda', lambda);

%% Look at factors
figure; SimpleWHPlot(W,H); title('SeqNMF reconstruction')
figure; SimpleWHPlot(W,H,X); title('SeqNMF factors, with raw data')

%% Procedure for choosing K
tic
Ws = {};
Hs = {};
numfits = 3; %number of fits to compare
for k = 1:10
    display(sprintf('running seqNMF with K = %i',k))
    for ii = 1:numfits
        [Ws{ii,k},Hs{ii,k}] = seqNMF(X,'K',k, 'L', L,'lambda', 0,'maxiter',30,'showplot',0);
        % note that max iter set low (30iter) for speed in demo (not recommended in practice)
    end
    inds = nchoosek(1:numfits,2);
    for i = 1:size(inds,1) % consider using parfor for larger numfits
            Diss(i,k) = helper.DISSX(Hs{inds(i,1),k},Ws{inds(i,1),k},Hs{inds(i,2),k},Ws{inds(i,2),k});
    end

end
%% Plot Diss and choose K with the minimum average diss.
figure,
plot(1:10,Diss,'ko'), hold on
h1 = plot(1:10,median(Diss,1),'k-','linewidth',2);
h2 = plot([3,3],[0,0.5],'r--');
legend([h1 h2], {'median Diss','true K'})
xlabel('K')
ylabel('Diss')

%% load example HVC calcium imaging data (from 6991FirstFewDaysForBatch)
clear all
display('Attempting to load MackeviciusData from seqNMF repository')
load MackeviciusData
display('loaded data')
%% break data into training set and test set
splitN = floor(size(NEURAL,2)*.75);
splitS = floor(size(SONG,2)*.75);
trainNEURAL = NEURAL(:,1:splitN);
trainSONG = SONG(:,1:splitS);
testNEURAL = NEURAL(:,(splitN+1):end);
testSONG = SONG(:,(splitS+1):end);
%% plot one example factorization
rng(235); % fixed rng seed for reproduceability
X = trainNEURAL;
K = 10;
L = 2/3; % units of seconds
Lneural = ceil(L*VIDEOfs);
Lsong = ceil(L*SONGfs);
shg
display('Running seqNMF on real neural data (from songbird HVC, recorded by Emily Mackevicius, Fee Lab)')
[W, H, ~,loadings,power]= seqNMF(X,'K',K,'L',Lneural,...
            'lambdaL1W', .1, 'lambda', .005, 'maxiter', 100, 'showPlot', 1,...
            'lambdaOrthoW', 0);
p = .05; % desired p value for factors

display('Testing significance of factors on held-out data')
[pvals,is_significant] = test_significance(testNEURAL,W,p);

W = W(:,is_significant,:);
H = H(is_significant,:);

% plot, sorting neurons by latency within each factor
[max_factor, L_sort, max_sort, hybrid] = helper.ClusterByFactor(W(:,:,:),1);
indSort = hybrid(:,3);
tstart = 180; % plot data starting at this timebin
figure; WHPlot(W(indSort,:,:),H(:,tstart:end), X(indSort,tstart:end), ...
    0,trainSONG(:,floor(tstart*SONGfs/VIDEOfs):end))
title('Significant seqNMF factors, with raw data')
figure; WHPlot(W(indSort,:,:),H(:,tstart:end), ...
    helper.reconstruct(W(indSort,:,:),H(:,tstart:end)),...
    0,trainSONG(:,floor(tstart*SONGfs/VIDEOfs):end))
title('SeqNMF reconstruction')

%% Procedure for choosing lambda
nLambdas = 20; % increase if you're patient
K = 10;
X = trainNEURAL;
lambdas = sort([logspace(-1,-5,nLambdas)], 'ascend');
loadings = [];
regularization = [];
cost = [];
for li = 1:length(lambdas)
    [N,T] = size(X);
    [W, H, ~,loadings(li,:),power]= seqNMF(X,'K',K,'L',Lneural,...
        'lambdaL1W', .1, 'lambda', lambdas(li), 'maxiter', 100, 'showPlot', 0);
    [cost(li),regularization(li),~] = helper.get_seqNMF_cost(X,W,H);
    display(['Testing lambda ' num2str(li) '/' num2str(length(lambdas))])
end
%% plot costs as a function of lambda
windowSize = 3;
b = (1/windowSize)*ones(1,windowSize);
a = 1;
Rs = filtfilt(b,a,regularization);
minRs = prctile(regularization,10); maxRs= prctile(regularization,90);
Rs = (Rs-minRs)/(maxRs-minRs);
R = (regularization-minRs)/(maxRs-minRs);
Cs = filtfilt(b,a,cost);
minCs =  prctile(cost,10); maxCs =  prctile(cost,90);
Cs = (Cs -minCs)/(maxCs-minCs);
C = (cost -minCs)/(maxCs-minCs);

clf; hold on
plot(lambdas,Rs, 'b')
plot(lambdas,Cs,'r')
scatter(lambdas, R, 'b', 'markerfacecolor', 'flat');
scatter(lambdas, C, 'r', 'markerfacecolor', 'flat');
xlabel('Lambda'); ylabel('Cost (au)')
set(legend('Correlation cost', 'Reconstruction cost'), 'Box', 'on')
set(gca, 'xscale', 'log', 'ytick', [], 'color', 'none')
set(gca,'color','none','tickdir','out','ticklength', [0.025, 0.025])

%% choose lambda=.005; run multiple times, see number of sig factors
loadings = [];
pvals = [];
is_significant = [];
X = trainNEURAL;
nIter = 20; % increase if patient
display('Running seqNMF multiple times for lambda=0.005')

for iteri = 1:nIter
    [W, H, ~,loadings(iteri,:),power]= seqNMF(X,'K',K,'L',Lneural,...
            'lambdaL1W', .1, 'lambda', .005, 'maxiter', 100, 'showPlot', 0);
    p = .05;
    [pvals(iteri,:),is_significant(iteri,:)] = test_significance(testNEURAL,W,p);
    W = W(:,is_significant(iteri,:)==1,:);
    H = H(is_significant(iteri,:)==1,:);
    [max_factor, L_sort, max_sort, hybrid] = helper.ClusterByFactor(W(:,:,:),1);
    indSort = hybrid(:,3);
    tstart = 300;
    clf; WHPlot(W(indSort,:,:),H(:,tstart:end), X(indSort,tstart:end), 0,trainSONG(:,floor(tstart*SONGfs/VIDEOfs):end))
    display(['seqNMF run ' num2str(iteri) '/' num2str(nIter)])
end
figure; hold on
h = histogram(sum(is_significant,2), 'edgecolor', 'w', 'facecolor', .7*[1 1 1]);
h.BinCounts = h.BinCounts/sum(h.BinCounts)*100;
xlim([0 10]);
xlabel('# significant factors')
ylabel('% seqNMF runs')

%% Plot factor-triggered song examples and rastors
addpath(genpath('misc_elm'));
figure; HTriggeredSpec(H,trainSONG,VIDEOfs,SONGfs,Lsong);
figure; HTriggeredRaster(H,trainNEURAL(indSort,:),Lneural);

%% Example parts-based and events-based factorizations
K = 3;
L = 50;
lambda =0;
X = NEURAL;

% run seqNMF with lambdaOrthoH -> events based
lambdaOrthoH = .1; % favor events-based (these can take any value, don't need to be zero and one)
lambdaOrthoW = 0;
display('Running seqNMF on simulated data, lambdaOrthoH -> events based')
figure;
[W,H] = seqNMF(X,'K',K, 'L', L,'lambda', lambda, ...
    'lambdaOrthoH', lambdaOrthoH, 'lambdaOrthoW', lambdaOrthoW);

% sort neurons and plot
[max_factor, L_sort, max_sort, hybrid] = helper.ClusterByFactor(W(:,:,:),1);
indSort = hybrid(:,3);
tstart = 180; % plot data starting at this timebin
WHPlot(W(indSort,:,:),H(:,tstart:end), X(indSort,tstart:end), ...
    1,trainSONG(:,floor(tstart*SONGfs/VIDEOfs):end)); title('lambdaOrthoH -> events based')

% run seqNMF with lambdaOrthoW -> parts based
figure;
lambdaOrthoH = 0;
lambdaOrthoW = 1; % favor parts-based (these can take any value, don't need to be zero and one)
display('Running seqNMF on simulated data, lambdaOrthoW -> parts based')
[W,H] = seqNMF(X,'K',K, 'L', L,'lambda', lambda, ...
    'lambdaOrthoH', lambdaOrthoH, 'lambdaOrthoW', lambdaOrthoW);

% sort neurons and plot
[max_factor, L_sort, max_sort, hybrid] = helper.ClusterByFactor(W(:,:,:),1);
indSort = hybrid(:,3);
WHPlot(W(indSort,:,:),H(:,:), X(indSort,:), ...
    1,trainSONG(:,:)); title('lambdaOrthoW -> parts based')

%% K sweep with masked cross-validation
nReps = 5; % increase if patient
Ks = 1:8; % increase if patient
L = 50;
X = NEURAL;
[N,T] = size(NEURAL);
RmseTrain = zeros(length(Ks), nReps);
RmseTest = zeros(length(Ks), nReps);
figure
[~,Kplot] = meshgrid(1:nReps, Ks);
Kplot = Kplot + rand(length(Ks), nReps)*.25-.125;
parfor K = Ks
    for repi = 1:nReps
        display(['Cross validation on masked test set; Testing K = ' num2str(K) ', rep ' num2str(repi)])
        rng('shuffle')
        M = rand(N,T)>.05; % create masking matrix (0's are test set, not used for fit)
        [W,H] = seqNMF(X,'K', K, 'L', L,'lambda', 0,'showPlot', 0, 'M', M);
        Xhat = helper.reconstruct(W,H);
        RmseTrain(K,repi) = sqrt(sum(M(:).*(X(:)-Xhat(:)).^2)./sum(M(:)));
        RmseTest(K,repi) = sqrt(sum((~M(:)).*(X(:)-Xhat(:)).^2)./sum(~M(:)));
    end
end

clf; scatter(Kplot(:), RmseTrain(:), 'r', 'markerfacecolor', 'flat');
hold on;
scatter(Kplot(:), RmseTest(:), 'b', 'markerfacecolor', 'flat');
plot(mean(RmseTrain,2), 'r')
plot(mean(RmseTest,2), 'b')
xlabel('K'); ylabel('RMSE')
legend('Train', 'Test', 'location', 'northwest')
drawnow; shg

%% Calculate the sequenciness score
% WARNING TAKES A WHILE
load MackeviciusData

nRepsShuff = 15;
nRepsColShuff = 15;  % just making an estimate, would need more to test sig
L = 20; % same as demo
K = 3;

X = NEURAL;
[N T] = size(X);

% do seqNMF
tmp = [];

parfor iteri = 1:nIter
    rng('shuffle')
    [~, ~, ~,~,tmp(iteri)] = seqNMF(X, 'L', L, 'K', K, 'lambda', 0, 'showPlot',0);
end

PEx = max(tmp);

% do seqNMF on shuffled data
PExShuff = [];
parfor repi = 1:nRepsShuff
    Xshuff = [];
    for ni = 1:N
        timeshuff = randperm(T);
        Xshuff(ni,:) = X(ni, timeshuff);
    end
    tmp = [];

    for iteri = 1:nIter
        rng('shuffle')
        [~, ~, ~,~,tmp(iteri)] = seqNMF(Xshuff, 'L', L, 'K', K, 'lambda', 0.0, 'showPlot',0);
    end
    PExShuff(repi) = max(tmp);
end

% do seqNMF on col shuffled data
PExColShuff = [];
parfor repi = 1:nRepsColShuff
    Xshuff = X(:,[1:L (L + randperm(T-L))]); % don't shuffle to first L bins... these cannot be explained by seqNMF
    tmp = [];

    for iteri = 1:nIter
        rng('shuffle')
        [~, ~, ~,~,tmp(iteri)] = seqNMF(Xshuff, 'L', L, 'K', K, 'lambda', 0.0, 'showPlot',0);
    end
    PExColShuff(repi) = max(tmp);
end
NoiseFloor = median(PExShuff);
SyncFloor = median(PExColShuff);
PAS = (PEx-SyncFloor)./...
    (PEx-NoiseFloor)

author: FeeLab

repo: seqNMF

save_path: github-repos/MATLAB/FeeLab-seqNMF


text:
function [ W, P, anim ]=generateW( anim, varargin )
% Generate projected data
%
% It uses the camera info and the 3D data to project to 2D
%
% USAGE
%  anim = anim.generateW()
%
% INPUTS
%  anim     - Animation object (help Animation for details)
%  varargin   - list of paramaters in quotes alternating with their values
%    P        - projection matrix. If not given, anim.P is used, if it does
%               not exist, it is computed from K,R,t
%    doFillP  - if set to true and anim.P is empty, it will fill anim.P
%    doFillW  - if set to true, will overwrite/create anim.
%
% OUTPUTS
%  W        - projected shape
%  P        - projection matrix
%  anim     - modified Animation
%
% EXAMPLE
%
% See also GENERATETOYANIMATION
%
% Vincent's Structure From Motion Toolbox      Version 3.1
% Copyright (C) 2008-2011 Vincent Rabaud.  [vrabaud-at-cs.ucsd.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the GPL [see external/gpl.txt]

[ P, doFillP, doFillW ] = getPrmDflt( varargin, ...
  { 'P' [] 'doFillP' false 'doFillW' false }, 1 );

W=[];

if isempty(anim.S)
  warning('S is not defined nor computable'); return;
end

if isempty(P)
  if isempty(anim.P); [ P anim ]=generateP(anim, doFillP);
  else P=anim.P;
  end
end

if ~isempty(anim.S)
  if size(P,2)==3
    % homography case
    W = bsxfun(@plus, multiTimes( P(:,1:2,:), anim.S, 1 ), P(:,3,:));
    W = reshape(normalizePoint(reshape(W,3,[]),3),2,anim.nPoint,...
      anim.nFrame);
  else
    % compute the projected points for normal camera and 3D points
    if size(anim.S,3)==1
      W = bsxfun( @plus, multiTimes( P(:,1:3,:), anim.S, 1 ), P(:,4,:) );
    else
      W = bsxfun( @plus, multiTimes( P(:,1:3,:), anim.S, 2 ), P(:,4,:) );
    end
    if anim.isProj; W = reshape(normalizePoint(reshape(W,3,[]),3),2,...
        anim.nPoint,anim.nFrame);
    else W = W(1:2,:,:);
    end
  end
end

% set missing values to NaN
if ~isempty(anim.mask); W(:,find(~anim.mask(:)))=NaN; end

if doFillW; anim.W=W; end
if doFillP; anim.P=P; end

author: vrabaud


text:
% x = blk2vec(X,nL)
%
% Converts a block diagonal matrix into a vector.
%
% **********  INTERNAL FUNCTION OF FROMPACK  **********

function x = blk2vec(X,nL)

%
% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA %

nblk = length(nL);
sumnk = 0;
x = [];
for k = 1:nblk
    nk = nL(k);
    x = [x; vec(X(sumnk+1:sumnk + nk,sumnk+1:sumnk + nk))] ;
    sumnk = sumnk + nk;
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% construct the Target Active Feature (TAF) network
% This net exploits the ridge loss and the hinge loss to calculate target
% active features.

% Input:
%  filter_sz    -   filter size [height width depth], target size in the
%                   input feature map
% Output:
%  IAF_net      -   the initialized IAF network

% By Xin Li, April 8 2018

function [TAF_net,filter_sz] = TAF_net_init(filter_sz)
rng('default');

channel=filter_sz(3);

rw=ceil(filter_sz(2)/2);
rh=ceil(filter_sz(1)/2);

fw=2*rw+1;
fh=2*rh+1;

filter_sz = [fh,fw,channel];

TAF_net=dagnn.DagNN();

%% conv layer 1-1 for regression
TAF_net.addLayer('conv11', dagnn.Conv('size', [fw,fh,channel,1],...
    'hasBias', true, 'pad',...
    [rh,rh,rw,rw], 'stride', [1,1]), 'input', 'conv11', {'conv11_f', 'conv11_b'});

f = TAF_net.getParamIndex('conv11_f') ;
TAF_net.params(f).value=single(randn(fh,fw,channel,1) /...
    sqrt(rh*rw*channel))/1e8;
TAF_net.params(f).learningRate=1;
TAF_net.params(f).weightDecay=1e3;

f = TAF_net.getParamIndex('conv11_b') ;
TAF_net.params(f).value=single(zeros(1,1));
TAF_net.params(f).learningRate=2;
TAF_net.params(f).weightDecay=1e3;

%%
TAF_net.addLayer('L2Loss',...
    RegressionL2Loss(),{'conv11','label_gaussian'},'objective_r');

% for the scale sensitive features, we first select the discriminative
% features with the regression loss then we select the scale sensitive from
% them with the ranking loss. The ranking loss is constructed in the file
% TAF_model.m

end

author: XinLi-zn

repo: TADT

save_path: github-repos/MATLAB/XinLi-zn-TADT


text:
classdef GeneralFormattingDirective < MBeautifier.Directive

    properties (SetAccess = immutable)
        Name = 'Format'
    end

    methods (Access = protected)
        function [success, toBeGarbadeCollected] = tryToupdateFromValue(obj, value)
            success = true;
            toBeGarbadeCollected = false;
            if strcmpi(value, 'on')
                obj.Values = {};
                toBeGarbadeCollected = true;
            elseif strcmpi(value, 'off')
                obj.Values = {'off'};
            else
                toBeGarbadeCollected = ~numel(obj.Values);
                success = false;
            end
        end
    end
end

author: davidvarga

repo: MBeautifier


text:
function [all_theta] = oneVsAll(X, y, num_labels, lambda)
%ONEVSALL trains multiple logistic regression classifiers and returns all
%the classifiers in a matrix all_theta, where the i-th row of all_theta
%corresponds to the classifier for label i
%   [all_theta] = ONEVSALL(X, y, num_labels, lambda) trains num_labels
%   logisitc regression classifiers and returns each of these classifiers
%   in a matrix all_theta, where the i-th row of all_theta corresponds
%   to the classifier for label i

% Some useful variables
m = size(X, 1);
n = size(X, 2);

% You need to return the following variables correctly
all_theta = zeros(num_labels, n + 1);

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: You should complete the following code to train num_labels
%               logistic regression classifiers with regularization
%               parameter lambda.
%
% Hint: theta(:) will return a column vector.
%
% Hint: You can use y == c to obtain a vector of 1's and 0's that tell use
%       whether the ground truth is true/false for this class.
%
% Note: For this assignment, we recommend using fmincg to optimize the cost
%       function. It is okay to use a for-loop (for c = 1:num_labels) to
%       loop over the different classes.
%
%       fmincg works similarly to fminunc, but is more efficient when we
%       are dealing with large number of parameters.
%
% Example Code for fmincg:
%
%     % Set Initial theta
%     initial_theta = zeros(n + 1, 1);
%
%     % Set options for fminunc
%     options = optimset('GradObj', 'on', 'MaxIter', 50);
%
%     % Run fmincg to obtain the optimal theta
%     % This function will return theta and the cost
%     [theta] = ...
%         fmincg (@(t)(lrCostFunction(t, X, (y == c), lambda)), ...
%                 initial_theta, options);
%

for c = 1:num_labels

% Set initial theta
initial_theta=zeros(n+1,1);
%initial_theta=zeros(num_labels, n + 1);

% Set options for fminunc
options = optimset('GradObj', 'on', 'MaxIter', 50);

% Run fmincg to obtain the optimal theta
% This function will return theta and the cost
[theta]=fmincg(@(t)(lrCostFunction(t,X, (y == c), lambda)), initial_theta, options);

% Updating theta to overall all_theta
if (c==1)
	history_theta=theta';
else
	history_theta=[history_theta; theta'];
end

end

all_theta=history_theta;

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

end

author: yhyap

repo: machine-learning-coursera


text:
function [ a, b ] = p39_lim ( )

%*****************************************************************************80
%
%% P39_LIM returns the integration limits for problem 39.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = 0.0;
  b = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function simpletrans
% 1D transport - modelling with extensions for decay and linear sorption
%    using mixing cell method
%
%   $Ekkehard Holzbecher  $Date: 2006/02/08 $
%--------------------------------------------------------------------------
T = 2;                     % maximum time [s]
L = 1;                     % length [m]
D = 0.1;                   % dispersivity [m*m/s]
v = 1;                     % velocity [m/s]
lambda = 1.2;              % decay constant [1/s]
R = 1;                     % retardation [1]
c0 = 0;                    % initial concentration [kg/m*m*m]
cin = 1;                   % inflow concentration [kg/m*m*m]

dtout = 0.05;              % output-timestep [s]
dxmax = 0.02;              % maximum grid spacing [m]
%------------------------ output parameters
gplot = 2;                 % =1: breakthrough curves; =2: profiles
gsurf = 0;                 % surface
gcont = 0;                 % =1: contours; =2: filled contours
ganim = 2;                 % animation

%------------------------ execution----------------------------------------

dtout = dtout/R;           % timestep reduction for retardation case
dx = dtout*v;              % grid spacing
K = 1;                     % K = reduction factor for grid spacing
if (dx>dxmax) K = ceil(dx/dxmax); end
dx = dx/K;                 % reduced grid spacing
dtadv=dtout/K;             % advection-timestep
N = ceil(L/dx);            % N = number of cells
x = linspace(0,(N-1)*dx,N);% nodes on x-axis
Neumann = D*dtadv/dx/dx;   % Neumann-number for dispersion
M = max (1,ceil(3*Neumann)); % M = reduction factor to fulfill Neumann-condition
Neumann = Neumann/M/R;     % reduced Neumann-number
dtdiff = dtadv/M;          % diffusion timestep
t = dtadv;

clear c c1 c2;
c(1:N) = c0; c1 = c;
k = 1; kanim = 1;
while (t < T/R)
    for i=1:M
        kinetics;          % decay (1. order kinetics)
        diffusion;         % diffusion
    end
    advection;             % advection
    if k >= K
        c = [c;c1]; k=0;
    end
    t = t + dtadv; k = k+1;
end
xlabel ('space'); ylabel ('concentration');

%-------------------- graphical output-------------------------------------

switch gplot
    case 1
        plot (c)        % breakthrough curves
        xlabel ('time'); ylabel ('concentration');
    case 2
        plot (x,c','--')% profiles
        xlabel ('space'); ylabel ('concentration');
end
if gsurf                % surface
    figure; surf (x,[0 t],c);
    xlabel ('space'); ylabel ('time'); zlabel('concentration');
end
if gcont figure; end
switch gcont
    case 1
        contour (c)     % contours
        grid on; xlabel ('space'); ylabel ('time');
    case 2
        contourf(c)     % filled contours
        colorbar; xlabel ('space'); ylabel ('time');
end
if (ganim)
    [FileName,PathName] = uiputfile('*.mpg');
    figure; if (ganim > 1) hold on; end
    for j = 1:size(c,1)
        axis manual;  plot (x,c(j,:),'r','LineWidth',2);
        YLim = [min(c0,cin) max(c0,cin)];
        legend (['t=' num2str(dtout*(j-1))]);
        Anim(j) = getframe;
        plot (x,c(j,:),'b','LineWidth',2);
    end
    mpgwrite (Anim,colormap,[PathName '/' FileName]);     % mgwrite not standard MATLAB
    movie (Anim,0);   % play animation
end

author: Sable

repo: mcbench-benchmarks


text:
function test03 ( dim_num, level_max )

%*****************************************************************************80
%
%% TEST03 call SPARSE_GRID_GL to create a Gauss-Legendre sparse grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer LEVEL_MAX, the level.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST03:\n' );
  fprintf ( 1, '  SPARSE_GRID_GL makes a sparse Gauss-Legendre grid.\n' );

  level_min = max ( 0, level_max + 1 - dim_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  LEVEL_MIN = %d\n', level_min );
  fprintf ( 1, '  LEVEL_MAX = %d\n', level_max );
  fprintf ( 1, '  Spatial dimension DIM_NUM = %d\n', dim_num );
%
%  Determine the number of points.
%
  point_num = sparse_grid_gl_size ( dim_num, level_max );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of unique points in the grid = %d\n', point_num );
%
%  Compute the weights and points.
%
  [ grid_weight, grid_point ] = sparse_grid_gl ( dim_num, level_max, point_num );
%
%  Print them out.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Grid weights:\n' );
  fprintf ( 1, '\n' );
  for point = 1 : point_num
    fprintf ( 1, '  %4d  %10f\n', point, grid_weight(point) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Grid points:\n' );
  fprintf ( 1, '\n' );
  for point = 1 : point_num
    fprintf ( 1, '  %4d', point );
    for dim = 1 : dim_num
      fprintf ( 1, '%10f', grid_point(dim,point) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
% JN Kather, extract microns per pixel MPP from Leica SVS
function currMPP = Info2MPP(currImInfo,channel)

    if isfield(currImInfo,'ImageDescription')
    currDescription = currImInfo(channel).ImageDescription;
    disp(currDescription);
    try
        st = strfind(currDescription,'MPP = ');         % find start
        currDescription = currDescription((st+6):end);  % strip before start
        en = strfind(currDescription,'|');      % find end
        if ~isempty(en)
        currDescription = currDescription(1:(en(1)-1));    % strip after end
        end
        currMPP = str2double(currDescription);  % get MPP (microns per pixel)
    catch
        warning('FAILED to parse description');
        currMPP = [];
    end
    else
        warning('no ImageDescription field');
        currMPP=[];
    end

end

author: jnkather

repo: MSIfromHE

save_path: github-repos/MATLAB/jnkather-MSIfromHE


text:
function[u,v,w]=sphere2uvw(lat,lon,dr,dth,dphi)
%SPHERE2UVW Converts a 3D spherical vector to a 3D Cartesian vector.
%
%   [U,V,W]=XYZ2SPHERE(LAT,LON,V1,V2,V3) converts a vector in spherical
%   coordinates with components V1, V2, and V3 located at point
%   (LAT, LON) into a Cartesian 3-vector with components U, V, and W.
%
%   LAT and LON are in degrees.
%
%   The vector in the spherical coordinate system has radial component
%   V1, longitudinal component V2, and latitudinal component V3.
%
%   The Cartesian vector [U,V,W] is in a reference frame with the
%   X-axis at zero degrees longitude and the Z-axis at the North Pole.
%
%   All input arguments should be arrays of the same size.
%
%   SPHERE2UVW is inverted by UVW2SPHERE.
%
%   See JSPHERE for related functions.
%
%   Usage: [x,y,z]=sphere2uvw(lat,lon,v1,v2,v3);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2007--2014 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(lat, '--t')
    sphere2uvw_test,return
end
[phi,theta]=jdeg2rad(lat,lon);

u      =  cos(phi).*cos(theta).*dr  -sin(theta).*dth -  sin(phi).*cos(theta).*dphi;
v      =  cos(phi).*sin(theta).*dr + cos(theta).*dth -  sin(phi).*sin(theta).*dphi ;
w      =  sin(phi).*dr                                + cos(phi).*dphi;

function[]=sphere2uvw_test

lat=[0  0  45         0        -90]';
lon=[0  90 0         90         0 ]';
v1= [1  1  1          0         0]';
v2= [0  0  0         -1         0]';
v3= [0  0  0          0         1]';
u=  [1  0  sqrt(2)/2  1         1 ]';
v=  [0  1  0          0         0]';
w=  [0  0  sqrt(2)/2  0         0]';

[u2,v2,w2]=sphere2uvw(lat,lon,v1,v2,v3);
tol=1e-6;
reporttest('SPHERE2XYZ example points',aresame(u,u2,tol)&&aresame(v,v2,tol)&&aresame(w,w2,tol))

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function Write_Energy(FileWritingSetup)
% Write_Energy(FileWritingSetup)
%
% FUNCTION:
%
%   This function takes the symbolic expressions for system energy,
%   position, and velocity, and uses them to write a function file.
%
%
% INPUTS: FileWritingSetup.
%
%   States = [2N x 2] cell array where the first column is the name of
%       the states, and the second column is comments about that states.
%       The index of a state in this array will match the index in the
%       dynamics function to be written. This function assumes a second
%       order system that is being written in first order form, so if
%       States{i,1} = z, then States{i+N,1} = dz; Furthermore, the output
%       (dStates) will have the form: dStates(i,:) = dz,
%       and dStates(i+N,:) = ddz.
%
%   Parameters = [K x 2] cell array with all of the input parameters. The
%       first column is for names, the second for comments
%
%   Energy = A struct with fields for the symbolic expressions for
%       position and velocity.
%
%   Directory = the name of the directory to write the files in.
%
%
% OUTPUTS:
%   This function will write the function kinematics()
%
%
% NOTE:
%
%   1)  It is assumed that the symbols in all of the inputs are self
%       consistent; that is - if Contacts{1,1} = H1, then there must be a field
%       in Dyn.(FileNames{i,1}) called H1.
%
%
% Written by Matthew Kelly
% December 6, 2013
% Cornell University
%
% See also DERIVE_EOM

States = FileWritingSetup.States;
Parameters = FileWritingSetup.Parameters;
Directory = FileWritingSetup.Directory;
Energy = FileWritingSetup.Energy;

CurrDir = cd;
if ~exist(Directory,'dir')
    mkdir(Directory);
end
cd(Directory);

fid = fopen('energy.m','w');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                         Function Header                             %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

fprintf(fid,'function Energy = energy(States, Parameters)\n');
fprintf(fid,'%% function Energy = energy(States, Parameters)\n');
fprintf(fid,'%%\n%% Computer Generated File -- DO NOT EDIT \n');
fprintf(fid,['%%\n%% This function was created by the function '...
    'Write_Energy()\n']);
fprintf(fid,['%% ' datestr(now) '\n%%\n']);
fprintf(fid,'%%  RETURNS: \n');
fprintf(fid,'%%   Energy = a struct with fields: \n');
fprintf(fid,'%%       Potential \n');
fprintf(fid,'%%       Kinetic \n');
fprintf(fid,'%%       Total \n');
fprintf(fid,'%% \n');
fprintf(fid,'%%\n');
fprintf(fid,'%% Matthew Kelly \n%% Cornell University \n%% \n\n');
fprintf(fid,'%% \n');
fprintf(fid,'%% See also DERIVE_EOM \n');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                           Parse Inputs                              %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

for idxState=1:size(States,1)
    fprintf(fid,[States{idxState,1} ' = States(:,' num2str(idxState) '); %% ' ...
        States{idxState,2} '\n']);
end %idxState
fprintf(fid,'\n');

for idxParam=1:size(Parameters,1)
    fprintf(fid,[Parameters{idxParam,1} ' = Parameters.' ...
        Parameters{idxParam,1} '; %% ' ...
        Parameters{idxParam,2} '\n']);
end %idxState
fprintf(fid,'\n');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                          Write Outputs                              %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

fprintf(fid,['Energy.Potential.m1 = ' ...
    vectorize(Energy.Potential.m1) ';\n']);
fprintf(fid,['Energy.Potential.m2 = ' ...
    vectorize(    Energy.Potential.m2) ';\n']);
fprintf(fid,['Energy.Potential.M = ' ...
    vectorize(Energy.Potential.M) ';\n']);
fprintf(fid,'Energy.Potential.Total = Energy.Potential.m1 + Energy.Potential.m2 + Energy.Potential.M;\n');
fprintf(fid,'\n');
fprintf(fid,['Energy.Kinetic.m1 = ' ...
    vectorize(Energy.Kinetic.m1) ';\n']);
fprintf(fid,['Energy.Kinetic.m2 = ' ...
    vectorize(Energy.Kinetic.m2) ';\n']);
fprintf(fid,['Energy.Kinetic.M = ' ...
    vectorize(Energy.Kinetic.M) ';\n']);
fprintf(fid,'Energy.Kinetic.Total = Energy.Kinetic.m1 + Energy.Kinetic.m2 + Energy.Kinetic.M;\n');
fprintf(fid,'\n');
fprintf(fid,'Energy.Total = Energy.Kinetic.Total + Energy.Potential.Total;\n');
fprintf(fid,'\n');
fprintf(fid,'end\n');
fclose(fid);

cd(CurrDir);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  ~~~~~~~~~~~~~~~~~~~ Sub Functions  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function characterExpression = vectorize(symbolicExpression)

        characterExpression = char(symbolicExpression);
        characterExpression = strrep(characterExpression,'*','.*');
        characterExpression = strrep(characterExpression,'/','./');
        characterExpression = strrep(characterExpression,'^','.^');
end

author: MatthewPeterKelly

repo: dscTutorials


text:
%CAVITYCONTROL demonstrates the control of lid-driven cavity flow using the
%Koopman-MPC framework described in
% "A data-driven Koopman model predictive control framework for nonlinear
% flows" H. Arbabi, M. Korda and I. Mezic

clc,clear
addpath('./thehood')
% quadratic programming solver
addpath('./thehood/qpOASES-3.1.0/interfaces/matlab')

if (exist('qpOASES_sequence','file') ~= 3)

    error(['You have to activate the MATLAB interface for qpOASES first:' ...
        ' unzip the qpOASES-3.1.0 in "thehood" folder, then '...
        'go to ".\thehood\qpOASES-3.1.0\interfaces\matlab" and run make.m,' ...
        'then run again'])

end

% simulation and flow parameters
SimPar.Reynolds = 13000;
SimPar.dt =.01;            % time step of the ODE solver
SimPar.N = 49;             % size of the computaional grid is (N+1)^2
SimPar.T = 0.2;            % length of each run - continuous-time cavity flow is approximated
                           % using T-time discrete map

load('CavityStateLibrary.mat','LimitCycle_Re13k','FixedPoint_Re10k','UnstableFixedPoint_Re13k')

% initial condition
 x0 = LimitCycle_Re13k(:,25);   % some point on the limit cyle

%% ************************PART I: Identification ********************** %%
%% STEP2: Identification via EDMD
% there are two options:
% 1- generate the data via runing the following program
% (for default values takes ~10hrs on a powerful desktop)
% DataFileName=GenerateCavityData;

% 2- download the data file cited in README.md and set
DataFileName = 'Cavity_data_4EDMD_0.mat';

% Generating Koopman-linear predictors with various number of measurements
CavitySystemID(DataFileName);

%  load any of the predictors
load('KoopmanLinSys_Re13_k50_random');    % predictor with k=50 random measurements

%% ************************PART II: control ********************** %%
% set up MPC controller
disp('setting up MPC controller ...')

n = size(A,1);  % state dimension
r = size(C,1);  % output dimension

% cost matrices
Qy = speye(r);  % output weight matrix
R = 0; % zero weight on input (there are input constraints)

% Prediction horizon
Tpred = 10;
Np = round(Tpred / SimPar.dt);

% State constraints
x_min = nan(n,1);
x_max = nan(n,1);

% Input constraints
umin = 11/13;
umax = 15/13;

% Precompute and possibly save MPC big matrices to speed construction of the controller.
% (This "lifting" has nothing to do with the koopman lifting)
[Ab, Bb] = createMPCmatrices(A,{B},Np);

% compute controller handle
tic
[~,~,mpcCont_lift]  = qpOases_MPC_controller(A,B,C,0,Qy,R,Qy,Np,umin, umax, x_min, x_max,'qpoases',[],[],[],[],[],Ab,Bb);
toc

% reference state
xref = FixedPoint_Re10k;    % unstable fixed point at Re=13k
xref_mpc = xref - x_mean; % mean subtracted!

% closed-loop simulation setup
% nonlinear solver
f = @(x,u)(NonlinearFlowSolver(x,CreateLidVelocity(u,SimPar.N),SimPar));

%% initialize
X = x0;
U = [];
cost = (CollectOutput(x0 - xref))'*Qy*(CollectOutput(x0 - xref));
Q_KE = getCostMatrix(SimPar.N); % the weight matrix to compute kinetic energy
KE_discrepancy = (x0 - xref)'*Q_KE*(x0 - xref);

% run to build the initial delay embeded state
for i = 1:nd

%     Nonlinear simulation with base input
    X = [ X, f(X(:,end),1) ];
    U = [U 1];
end

% now build the initial embedded state
disp('initialization complete')

%% closed-loop simulation

Tsim = 100;
Nsim = Tsim / SimPar.T;

for i = 1:Nsim
    if mod(i,50)==0
        fprintf('Closed-loop simulation, %f %% completed \n', 100 * i / Nsim)
    end
    % create reference output
    yref = CollectOutput(xref_mpc(:,end));

    % build the state of the Koopman predictor
    y = CollectOutput(  bsxfun(@minus,X(:,end-nd+1:end),x_mean) ); % collect the measurements
    y = DelayEmbed(y,nd);   % delay embed the measurements
    ue = DelayEmbed(U(:,end-nd+1:end-1),nd-1);  % delay embed the input
    znow = [y;ue;KE_embed(y);1];  % form the state (g(zeta) in the paper)

    % compute control
    u = mpcCont_lift(znow,yref);

    % Nonlinear simulation
    X = [ X, f(X(:,end),u(:,1)) ];
    U = [U u(:,1)];

    % compute the tracking error and kinetic energy of state discrepancy
    ynow = CollectOutput(  bsxfun(@minus,X(:,end),x_mean) );
    cost = [ cost, (ynow-yref)'*Qy*(ynow-yref) ];
    KE_discrepancy = [KE_discrepancy, (X(:,end) - xref)'*Q_KE*(X(:,end) - xref)];

end

%% plots
set(0,'defaultTextInterpreter','latex', ...
    'defaultLegendInterpreter','latex', ...
    'defaultAxesTickLabelInterpreter','latex');
t = (0:Nsim)*SimPar.T;
figure(20),clf
subplot(2,2,1)
plot(t,U(1,nd:end),'linewidth',2);
xlabel('$t$','fontsize',12);
title('input')
ylim([umin umax])
set(gca,'YTick',linspace(umin,umax,5))

subplot(2,2,2)
plot(t,KE_discrepancy); hold on
title('kinetic energy of discrepency')

subplot(2,2,3)
plot(t,cost); hold on
title('control tracking error')

subplot(2,2,4)
PlotVorticity(X(:,end)-xref);
hold on
[ Grid ] = CavityGridOperators( SimPar.N );
plot(CollectOutput(Grid.xx),CollectOutput(Grid.yy),'x')
axis square
axis([-1 1 -1 1])
title('vorticity discrepancy at final time and sensor locations')

author: arbabiha


text:
function model = addCOBRAVariables(model, idList, varargin)
% Add a Variable to the COBRA model provided. These variables can be
% referred to by additional constraints added via addCOBRAConstraint. Their
% IDs have to be mutually exclusive with the IDs in model.rxns
% USAGE:
%    model = addCOBRAVariables(model, idList, varargin)
%
% INPUTS:
%    model:         model structure
%    idList:        cell array of ids
%
%    varargin:
%                   * lb:               The lower bounds of the variables
%                                       Has to be of the same size as
%                                       idlist
%                                       (default: -1000)
%                   * ub:               The upper bounds of the variables.
%                                       Has to be of the same size as
%                                       idlist
%                                       (default: 1000)
%                   * c:                The objective coefficient of the variables
%                                       Has to be of the same size as idlist
%                                       (default: 0)
%                   * Names:            Descriptive names of the variables.
%                                       (default: idList)
% OUTPUT:
%    model:         constrained model containing the new variables in the
%                   respective fields (`E`,`evars`,`evarub`,`evarlb`,`evarc`)
%                    * `.E` - The coefficient matrix for coefficients of metabolites for the variables
%                    * `.evars` - The variable IDs
%                    * `.evarlb` - The variable lower bounds
%                    * `.evarub` - The variable upper bounds
%                    * `.evarc` - The variable objective coefficient
%                    *  Optional: `D`, the matrix conatining coefficients for additional variables.
% NOTE:
%    This function will, if not present create the `E`, `evars`, `evarlb`, `evarub`, `evarc` and,
%    if additional Constraints (i.e. the `C` matrix) is present, the `D`
%    field as defined in the model field definitions. The `E` matrix containts
%    the Variable coefficients referring to metabolites, while the `D`
%    field contains the variable coefficients referring to additional
%    Constraints from `C`. `evars` contains the IDs of the additional
%    variables, with one entrey per column of `E`. `evarlb` and `evarub`
%    represent the lower and upper bounds of the variables respectively.
%    `evarc` are objective coefficients of the variables.
%
% Author: Thomas Pfau, Oct 2018

if ischar(idList)
    %This is an individual element.
    idList = {idList};
end

p = inputParser();
p.CaseSensitive = false;
p.addRequired('model',@isstruct);
p.addRequired('idList',@iscell);
p.addParameter('lb',columnVector(-1000*ones(size(idList))),@(x) isnumeric(x) && (numel(x) == numel(idList) || numel(x) == 1));
p.addParameter('ub',columnVector(1000*ones(size(idList))),@(x) isnumeric(x) && (numel(x) == numel(idList) || numel(x) == 1));
p.addParameter('c',columnVector(zeros(size(idList))),@(x) isnumeric(x) &&( numel(x) == numel(idList) || numel(x) == 1));
p.addParameter('Names',idList, @(x) (iscell(x) && numel(x) == numel(idList)) || ischar(x) && numel(idList == 1));
p.parse(model,idList,varargin{:});

lb = p.Results.lb;
if numel(lb) == 1
    lb = ones(numel(idList),1)*lb;
end

ub = p.Results.ub;
if numel(ub) == 1
    ub = ones(numel(idList),1)*ub;
end

c = p.Results.c;
if numel(c) == 1
    c = ones(numel(idList),1)*c;
end

if any(ismember(model.rxns,idList))
    duplicateIDs = unique(model.rxns(ismember(model.rxns,idList)));
    error('The following IDs are already IDs of reactions:\n%s\n',strjoin(duplicateIDs,'\n'));
end

% no variables exist yet, so we need to build the fields.
if ~isfield(model,'evars')
    % we need to create new fields
    if ~isfield(model,'ctrs')
        fieldsToCreate = {'E','evars','evarlb','evarub','evarc'};
    else
        fieldsToCreate = {'E','evars','evarlb','evarub','evarc','D'};
    end
    model = createEmptyFields(model,fieldsToCreate);
else
    % check, that the IDs are valid.
    if any(ismember(model.evars,idList))
        duplicateIDs = unique(model.evars(ismember(model.evars,idList)));
        error('The following IDs are already IDs of variables:\n%s\n',strjoin(duplicateIDs,'\n'));
    end
end

% this is all checks done, now lets get to work.
[~,nVars] = size(model.E);

if ~isfield(model,'evarNames') && ~any(strcmp(p.UsingDefaults,'Names'))
    % if the field does not yet exist, and we got actual names, we need to
    % create the field.
    model = createEmptyFields(model,{'evarNames'});
end

% extend the IDs.
model.evars = [model.evars; columnVector(idList)];
model.evarlb = [model.evarlb; columnVector(lb)];
model.evarub = [model.evarub; columnVector(ub)];
model.evarc = [model.evarc; columnVector(c)];

if isfield(model,'evarNames')
    model.evarNames = [model.evarNames; columnVector(p.Results.Names)];
end

% extend other fields as necessary.
model = extendModelFieldsForType(model,'evars','originalSize',nVars);

author: opencobra

repo: cobratoolbox


text:
% NCGEOVARIABLE Provide advanced access to variables and their related
% dimensions.
%
% NCGEOVARIABLE is used to retrieve data for a given variable as well as the
% variables associated coordinate dimensions.
%
% Methods:
%
% For more information on the methods use 'help' or 'doc'.  For example:
%   >> help ncgeovariable
%   >> doc ncgeovariable
%
% Example of use:
%  ds = ncgeodataset('http://dods.mbari.org/cgi-bin/nph-nc/data/ssdsdata/deployments/m1/200810/OS_M1_20081008_TS.nc');
%  vn = ds.variable('TEMP');    % returns an ncvariable
%  vg = ds.geovariable('TEMP'); % returns an ncgeovariable
%  tn = vn.data([1 1 1 1], [100 5 1 1]);
%  tg = vg.data(1:100,1:5,1,1);
%  % Look at properties
%  vg.name
%  vg.axes
%  vg.getaxesorder
%  vg.extent
%  datestr(vg.timeextent)
%  vg.grid_interop
%  tw = vg.timewindow([2008 10 10 0 0 0], [2008 20 11 0 0 0]);
%  datestr(tw.time)
%  vg.getlonname
%  vg.getlatname
%
%  vg.getxname
%  vg.getyname
%
% NCTOOLBOX (https://github.com/nctoolbox/nctoolbox
% https://github.com/nctoolbox/nctoolbox/wiki/ncgeovariable)
% See also ncvariable, ncgeodataset, ncgeodataset.geovariable
classdef ncgeovariable < ncvariable

    properties (SetAccess = private)
        %         dataset          % ncdataset instance.
    end

    properties (Dependent = true)
        %         name            % The string variable name that this object represents.
        %         attributes      % The attributes associated with this object.
        %         axes            % The coordiante vaaibles associated with this object.
    end

    properties (SetAccess = private, GetAccess = protected)
        %         variable        % ucar.nc2.Variable instance. Represents the data
        %         axesVariables    % ucar.nc2.Variable instance. Represents the data.

    end

    properties (SetAccess = private, GetAccess = private)
        %         variable        % ucar.nc2.Variable instance. Represents the data
        %         axesVariables    % ucar.nc2.Variable instance. Represents the data.
        axes_info % list of axes names and the dimensions in one cell dict
    end

    methods

        %%
        function obj = ncgeovariable(src, variableName, axesVariableNames)
            % NCGEOVARIABLE.NCGEOVARIABLE  Constructor.
            %
            % Use as:
            %    v = ncgeovariable(src, variableName, axesVariableNames)
            %
            % Preffered Usage:
            %    geo = ncgeodataset(...);
            %    v = geo.geovariable(variableName); % uses ncgeodataset.geovariable
            % Alternate Usage:
            %    geo = ncgeodataset(...);
            %    v = ncgeovariable(geo,variableName,geo.axes(variableName))
            obj = obj@ncvariable(src, variableName, axesVariableNames);

        end % ncgeovariable end

        %%
        function a = get.axes_info(src)
            switch length(src.size)
                case 1
                    a = fieldnames( src.grid_interop(1,1,1));
                    b = fieldnames( src.grid(1,1,1));
                case 2
                    a = fieldnames( src.grid_interop([1,1],[1,1],[1,1]));
                    b = fieldnames( src.grid([1,1],[1,1],[1,1]));
                case 3
                    a = fieldnames( src.grid_interop([1,1,1], [1,1,1], [1,1,1]));
                    b = fieldnames( src.grid([1,1,1], [1,1,1], [1,1,1]));
                case 4
                    a = fieldnames( src.grid_interop([1,1,1,1], [1,1,1,1], [1,1,1,1]));
                    b = fieldnames( src.grid([1,1,1,1], [1,1,1,1], [1,1,1,1]));
                case 5
                    a = fieldnames( src.grid_interop([1,1,1,1,1], [1,1,1,1,1], [1,1,1,1,1]));
                    b = fieldnames( src.grid([1,1,1,1,1], [1,1,1,1,1], [1,1,1,1,1]));
                case 6
                    a = fieldnames( src.grid_interop([1,1,1,1,1,1], [1,1,1,1,1,1], [1,1,1,1,1,1]));
                    b = fieldnames( src.grid([1,1,1,1,1,1], [1,1,1,1,1,1], [1,1,1,1,1,1]));
            end
            len = length(a);
            for i = 1:len
                a{i, 2} = src.dataset.size(b{i});
            end

        end

        function e = top_bot(src)
            % NCGEOVARIABLE.top_bot - Find top and bottom layers
            % Usage e = var.top_bot;
            % returns [top_index bottom_index]
            s = src.size;
            lens = length(s);
            first = ones(1,lens);
            if lens == 3
                s(2) = 1;
                s(3) = 1;
            elseif lens==4
                s(1) = 1;
                s(3) = 1;
                s(4) = 1;
            else
            end
            stride = first;
            %             switch length(s) % hopefully this speeds up the grid_interop call when time is involved
            %                 case 1
            %                     g = src.grid_interop(:);
            %                 case 2
            %                     g = src.grid_interop(:,:);
            %                 case 3
            %                     g = src.grid_interop(1,:,:);
            %                 case 4
            g = src.grid_interop(first, s, stride);
            %             end
            if g.z(1) <= g.z(2)
                e=[length(g.z) 1];
            else
                e=[1 length(g.z)];
            end
        end

        function e = extent(src)
            % NCGEOVARIABLE.extent - Function to find geographic bounding box coordinates.
            % Usage e = var.extent;
            s = src.size;
            lens = length(s);
            first = ones(1,lens);
            if lens > 1
                s(1) = 1;
                if lens > 3;
                    s(2) = 1;
                end
            end
            stride = first;
            %             switch length(s) % hopefully this speeds up the grid_interop call when time is involved
            %                 case 1
            %                     g = src.grid_interop(:);
            %                 case 2
            %                     g = src.grid_interop(:,:);
            %                 case 3
            %                     g = src.grid_interop(1,:,:);
            %                 case 4
            g = src.grid_interop(first, s, stride);
            %             end
            e.lon = [min(min(g.lon)) max(max(g.lon))];
            e.lat = [min(min(g.lat)) max(max(g.lat))];
        end % end extent

        function te = timeextent(src)
            % NCGEOVARIABLE.timeextent - Function to find the start and stop time of the variable.
            % Usage e = var.timeextent;
            s = src.size;
            lens = length(s);
            first = ones(1, lens);
            stride = first;
            last = first;
            last(1) = s(1);

            %                     g1 = src.grid_interop(first, first, stride);
            %                     g2 = src.grid_interop(s, s, stride);
            g = src.grid_interop(first, last, stride);

            %             te = [g1.time g2.time];
            if isfield(g, 'time')
                te = [min(g.time) max(g.time)];
            else
                error('NCGEOVARIABLE:TIMEEXTENT',...
                    'There appears to be no time axis associated with the variable.');
            end

        end % end timeextent

        function ig = grid_interop(src, first, last, stride, ignore)
            % NCGEOVARIABLE.GRID_INTEROP - Method to get the coordinate variables and their data as a
            % a structure with standardized field names for lat, lon, time, and z. Other coordinate variables
            % that are not recognized by netcdf-java as the previous four types have field names directly
            % taken from their variable names.
            % Usage: >> gridstruct = geovar.grid_interop(1,:,:,1:2:50);
            if nargin < 5
                ignore = {};
            end
            g = src.grid(first, last, stride);
            names = fieldnames(g);

            for i = 1:length(names); % loop through fields returned by grid
                tempname = names{i};
                javaaxisvar  =   src.dataset.netcdf.findVariable(tempname);
                try
                    % Expecting a uvar.nc2.dataset.CoordinateAxis here
                    type = char(javaaxisvar.getAxisType());
                catch me
                    % If it's not a CoordinateAxis set type to an empty
                    % string
                    type = '';
                end
                if isempty(type)
                    ig.(tempname) = g.(tempname);
                else
                    switch type
                        case 'Height'
                            if ~ismember('z', ignore)
                                pos_z = char(javaaxisvar.getPositive());
                                if strcmpi(pos_z, 'POSITIVE_DOWN')
                                    tmp = g.(tempname);
                                    ig.z = tmp.*-1; %adjust for positive direction
                                elseif strcmpi(pos_z, 'down')
                                    tmp = g.(tempname);
                                    ig.z = tmp.*-1; %adjust for positive direction
                                else
                                    ig.z = g.(tempname);
                                end
                            end
                        case 'GeoZ'
                            if ~ismember('z', ignore)
                                pos_z = char(javaaxisvar.getPositive());
                                z_sn = src.dataset.attribute(tempname, 'standard_name');
                                k = strfind(z_sn, 'ocean_s');
                                switch isempty(k)
                                    case 0
                                        %                       n = max(size(src.size));
                                        %                       trange = java.util.ArrayList(n);
                                        %                       zrange = trange;
                                        %                       xrange = trange;
                                        %                       yrange = trange;
                                        trange = ucar.ma2.Range(first(1) - 1, last(1) - 1, stride(1));
                                        zrange = ucar.ma2.Range(first(2) - 1, last(2) - 1, stride(2));
                                        xrange = ucar.ma2.Range(first(4) - 1, last(4) - 1, stride(4));
                                        yrange = ucar.ma2.Range(first(3) - 1, last(3) - 1, stride(3));
                                        %                       coordinates.GridCoordSys.getVerticalTransform.getCoordinateArray
                                        %                          temp = src.variable.getCoordinateSystems();
                                        griddataset = ucar.nc2.dt.grid.GridDataset.open(src.dataset.location);
                                        grid = griddataset.findGridByName(src.name);
                                            if ( length(grid) ~= 0 )
                                            grid = grid.getCoordinateSystem();
                                            subgrid = grid.getVerticalTransform();
                                        else
                                            if ~ismember(tempname, ignore)
                                                ig.(tempname) = g.(tempname);
                                            end
                                            continue % skip to next var
                                        end
                                        subgrid = subgrid.subset(trange, zrange, yrange, xrange); %tempfortesting arc 10/18
                                        %
                                        % It looks like this dataset (http://geoport.whoi.edu/thredds/dodsC/usgs/vault0/models/examples/bora_feb.nc)
                                        % works when the vertical transform isnt subset. NJTBX uses the same methodology but subsets the resulting
                                        % z coordinate field in matlab after calling the getCoordinateArray method, not before in the java.
                                        % Should I implement a try and then default to subsetting afterwards on catch? Or should I just always do it
                                        % after? I am choosing the former. Anyone have any thoughts on the matter? -acrosby
                                        try
                                            try
                                                try
                                                    for q = first(1):stride(1):last(1)
                                                        %                                                     grid = griddataset.findGridByName(src.name);
                                                        %                                                     subgrid = grid.getCoordinateSystem();
                                                        array = subgrid.getCoordinateArray(q-1);
                                                        ig.z(q, :, :, :) = array.copyToNDJavaArray();
                                                    end
                                                catch me
                                                    c = 1;
                                                    for q = first(1):stride(1):last(1)
                                                        %                                                     grid = griddataset.findGridByName(src.name);
                                                        %                                                     grid = grid.getCoordinateSystem();
                                                        subgrid = grid.getVerticalTransform();

                                                        array = subgrid.getCoordinateArray(q-1); % Issue 27 is failing here...
                                                        z(c, :, :, :) = array.copyToNDJavaArray();
                                                        c = c + 1;
                                                    end
                                                    ig.z = z(:, first(2):stride(2):last(2),  first(3):stride(3):last(3),  first(4):stride(4):last(4));
                                                end
                                            catch me
                                                array = subgrid.getCoordinateArray(0);
                                                ig.z = array.copyToNDJavaArray();
                                            end
                                        catch me
                                            disp('Could you please add the code you are trying to run to Issue 27 at the nctoolbox issue tracking site.');
                                            web http://code.google.com/p/nctoolbox/issues/detail?id=27
                                            me.throw()
                                            % me.error('There is a problem applying the vertical coordinate tranform and subsetting the resuting values.');
                                        end

                                    otherwise

                                        if strcmp(pos_z, 'POSITIVE_DOWN')
                                            tmp = g.(tempname);
                                            ig.z = tmp.*-1; %adjust for positive direction
                                        else
                                            ig.z = g.(tempname);
                                        end
                                end

                            end
                        case 'Time'
                            if ~ismember('time', ignore)
                                tmp = g.(tempname);
                                t_converted = src.dataset.time(tempname, tmp);
                                ig.time = t_converted;
                            end
                            %                     case 'RunTime'
                            %                       tmp = obj.dataset.data(name, vFirst, vLast, vStride);
                            %                       t_converted = obj.dataset.time(name, tmp);
                            %                       data.(type) = t_converted;

                        case 'Lon'
                            if ~ismember('lon', ignore)
                                tmp = g.(tempname);
                                if max(max(tmp)) > 360
                                    if min(min(tmp)) > 0
                                        tmp(tmp>360) = tmp(tmp>360)-360;
                                        tmp(tmp>360) = tmp(tmp>360)-360;
                                        tmp(tmp>360) = tmp(tmp>360)-360;
                                        tmp(tmp>180) = tmp(tmp>180)-360;
                                    else
                                        error('NCGEOVARIABLE:GEOIJ',...
                                            'Longitude contains values that follow both -180/180 and 0/360+ conventions; can not subset.');
                                    end
                                elseif min(tmp) <~ 0
                                    % Do nothing, we assume that input is in -180/180 conventions. And this means that the data is too.
                                elseif max(tmp) > 180
                                    tmp = tmp - 360;
                                end
                                ig.lon = double(tmp);
                            end

                        case 'Lat'
                            if ~ismember('lat', ignore)
                                ig.lat = double(g.(tempname));
                            end

                        case 'GeoY'
                            %ig.y = g.(tempname);
                            if exist('griddataset', 'var')
                                grid = griddataset.findGridByName(src.name);
                                grid = grid.getCoordinateSystem();
                            else
                                griddataset = ucar.nc2.dt.grid.GridDataset.open(src.dataset.location);
                                grid = griddataset.findGridByName(src.name);
                                grid = grid.getCoordinateSystem();
                            end
                            try
                                %ig.y = grid.getYHorizAxis;
                                %ig.x = grid.getXHorizAxis;
                                [x, y] = meshgrid(g.x,g.y);
                                s = size(x);
                                x = reshape(x, [1 numel(x)]);
                                y = reshape(y, [1 numel(y)]);
                                tempXY = [x; y];
                                projection = grid.getProjection();
                                tempLatLon = projection.projToLatLon(tempXY);
                                if ~isfield(ig, 'lat')
                                    if ~ismember('lat', ignore)
                                        ig.lat = reshape(tempLatLon(1,:), s);
                                    end
                                end
                                if ~isfield(ig, 'lon')
                                    if ~ismember('lon', ignore)
                                        ig.lon = reshape(tempLatLon(2,:), s);
                                    end
                                end
                            catch me
                            end

                        case 'GeoX'
                            %                             ig.x = g.(tempname);
                            %                             if exist('griddataset', 'var')
                            %                                 grid = griddataset.findGridByName(src.name);
                            %                                 grid = grid.getCoordinateSystem();
                            %                             else
                            %                                 griddataset = ucar.nc2.dt.grid.GridDataset.open(src.dataset.location);
                            %                                 grid = griddataset.findGridByName(src.name);
                            %                                 grid = grid.getCoordinateSystem();
                            %
                            %                             end
                            %                             try
                            %                                 ig.x = grid.getXHorizAxis;
                            %                                 projection = grid.getProjection();
                            %                                 ig.lon = projection.projToLatLon();
                            %                             catch me

                        otherwise
                            if ~ismember(tempname, ignore)
                                ig.(tempname) = g.(tempname);
                            end
                    end % end switch on type
                end % end is type empty or not if statement
            end % end loop through field names

        end % grid_interop end

        function tw = timewindow(src, varargin)
            % NCGEOVARIABLE.TIMEWINDOW - Function to pull the time coordinates within the specified
            % start and stop times from the variable object.
            % Usage: >> time = geovar.timewindow([2004 1 1 0 0 0], [2005 12 31 0 0 0]);
            %              >> time = geovar.timewindow(731947, 732677);
            if nargin < 3
                d = src.timewindowij(varargin{1});
            elseif nargin == 3
                d = src.timewindowij(varargin{1}, varargin{2});
            else
                error('NCGEOVARIABLE:TIMEWINDOW',...
                    'Too many input arugments.');
            end
            tw = d;
        end

        function tv = gettimevar(src)
            % NCGEOVARIABLE.gettimevar() gets the time variable associated with the ncgeovariable.
            tn = src.gettimename();
            tv = src.dataset.geovariable(tn);
        end

        function [lv, type] = getlonvar(src)
            % NCGEOVARIABLE.getlonvar() gets the longitude variable assiciated with the ncgeovariable.
            tn = src.getlonname();
            type = 'Lon';
            if isempty(tn)
                tn = src.getxname();
                type = 'GeoX';
            end
            lv = src.dataset.geovariable(tn);
        end

        function [lv, type] = getlatvar(src)
            % NCGEOVARIABLE.getlatvar() gets the latitude variable associated with the ncgeovariable.
            type = 'Lat';
            tn = src.getlatname();
            if isempty(tn)
                tn = src.getyname();
                type = 'GeoY';
            end
            lv = src.dataset.geovariable(tn);
        end

        function tn = gettimename(src)
            % NCGEOVARIABLE.gettimename() gets the name if the time variable associated with the ncgeovariable
            for i = 1:length(src.axes)
                tempname = src.axes{i};
                javaaxisvar  =   src.dataset.netcdf.findVariable(tempname);
                type{i} = char(javaaxisvar.getAxisType());
            end
            match = strcmp('Time', type);
            tn = src.axes(match);
        end

        function ln = getlonname(src)
            % NCGEOVARIABLE.getlonname() gets the name of the longitude variable associated with the ncgeovariable.
            for i = 1:length(src.axes)
                tempname = src.axes{i};
                javaaxisvar  =   src.dataset.netcdf.findVariable(tempname);
                type{i} = char(javaaxisvar.getAxisType());
            end
            match = strcmp('Lon', type);
            ln = src.axes(match);
        end

        function ln = getlatname(src)
            % NCGEOVARIABLE.gelatname() gets the naem of the latitude variable associated with the ncgeovariable
            for i = 1:length(src.axes)
                tempname = src.axes{i};
                javaaxisvar  =   src.dataset.netcdf.findVariable(tempname);
                type{i} = char(javaaxisvar.getAxisType());
            end
            match = strcmp('Lat', type);
            ln = src.axes(match);
        end

        function ln = getxname(src)
            % NCGEOVARIABLE.getxname() gets the name of the x variable associated with the ncgeovariable.
            for i = 1:length(src.axes)
                tempname = src.axes{i};
                javaaxisvar  =   src.dataset.netcdf.findVariable(tempname);
                type{i} = char(javaaxisvar.getAxisType());
            end
            match = strcmp('GeoX', type);
            ln = src.axes(match);
        end

        function ln = getyname(src)
            % NCGEOVARIABLE.getyname() gets the name of the y variable assoxiated with the ncgeovariable
            for i = 1:length(src.axes)
                tempname = src.axes{i};
                javaaxisvar  =   src.dataset.netcdf.findVariable(tempname);
                type{i} = char(javaaxisvar.getAxisType());
            end
            match = strcmp('GeoY', type);
            ln = src.axes(match);
        end

        function s = gettimedata(src, start, last, stride)
            % NCGEOVARIABLE.gettimedata()
            %             var = src.gettimevar;
            %             tn = var.data(start, last, stride);
            %             tn = var.dataset.time(src.gettimename, tn);

            v = src.gettimevar;
            sz = src.size();
            timesize = v.size();
            timelocation = find(sz==timesize(1));

            timestart = start(timelocation);
            timelast = last(timelocation);
            timestride = stride(timelocation);

            s = v.data(timestart:timestride:timelast);

            s = src.dataset.time(v.name, s);
        end

        function [s, type] = getlondata(src, start, last, stride)
            % NCGEOVARIABLE.getlondata() gets the longitude data associated with the ncgeovariable.
            [v, type] = src.getlonvar;
            sz = src.size();
            lonsize = v.size();
            lonlocation = find(sz==lonsize(1));
            if length(lonsize) == 1
                lonstart = start(lonlocation);
                lonlast = last(lonlocation);
                lonstride = stride(lonlocation);
            else
                lonstart = start(lonlocation:lonlocation+1);
                lonlast = last(lonlocation:lonlocation+1);
                lonstride = stride(lonlocation:lonlocation+1);
            end
            switch length(lonsize)
                case 1
                    s = v.data(lonstart:lonstride:lonlast);
                case 2
                    s = v.data(lonstart(1):lonstride(1):lonlast(1),lonstart(2):lonstride(2):lonlast(2));
                case 3
                    s = v.data(lonstart(1):lonstride(1):lonlast(1), lonstart(2):lonstride(2):lonlast(2), lonstart(3):lonstride(3):lonlast(3));
                case 4
                    s = v.data(lonstart(1):lonstride(1):lonlast(1), lonstart(2):lonstride(2):lonlast(2), lonstart(3):lonstride(3):lonlast(3),...
                        lonstart(4):lonstride(4):lonlast(4));
            end
            %               s = v.data(lonstart, lonlast, lonstride);
        end

        function [s, type] = getlatdata(src, start, last, stride)
            % NCGEOVARIABLE.gelatdata() gets the latitude data associated with the ncgeovariable.
            [v, type] = src.getlatvar;
            sz = src.size();
            latsize = v.size();
            latlocation = find(sz==latsize(1));
            if length(latsize) == 1
                latstart = start(latlocation);
                latlast = last(latlocation);
                latstride = stride(latlocation);
            else
                latstart = start(latlocation:latlocation+1);
                latlast = last(latlocation:latlocation+1);
                latstride = stride(latlocation:latlocation+1);
            end
            switch length(latsize)
                case 1
                    s = v.data(latstart:latstride:latlast);
                case 2
                    s = v.data(latstart(1):latstride(1):latlast(1),latstart(2):latstride(2):latlast(2));
                case 3
                    s = v.data(latstart(1):latstride(1):latlast(1), latstart(2):latstride(2):latlast(2), latstart(3):latstride(3):latlast(3));
                case 4
                    s = v.data(latstart(1):latstride(1):latlast(1), latstart(2):latstride(2):latlast(2), latstart(3):latstride(3):latlast(3),...
                        latstart(4):latstride(4):latlast(4));
            end
            %                s = v.data(latstart, latlast, latstride);
        end

        function vec = getvectors(src, u, v, alpha)
            % NCGEOVARIABLE.getvectors - Function to average and rotate
            % u and v ROMS vectors on to a common grid (grid is based on
            % the coordinates of the geovariable that "getvectors" is called
            % from.)
            % Use:
            %     >> uvar = nc.geovariable('u');
            %     >> vvar = nc.geovariable('v');
            %     >> hvar = nc.geovariable('depth');
            %     >> avar = nc.geovariable('angle');
            %     >> U = hvar.getvectors(uvar, vvar, avar); % complex matrix is returned
            if nargin < 4
                vec = complex_ongrid(src, u, v);
            elseif nargin == 4
                vec = complex_ongrid(src, u, v, alpha);
            else
                error(['NCTOOLBOX:ncgeovariable:getvectors'], ...
                    'Incorrect number of arguments for getvectrs()');
            end
        end

        %% These functions would rather output multiple outputs instead of struct, must reconcile
        %     with the subsref in either ncgeovariable or ncvariable. Wait, why, then, does geoij work???
        function d = timewindowij(src, varargin)
            % NCGEOVARIABLE.TIMEWINDOWIJ - Function to get indices from start and stop times for sub-
            % setting. TODO: There must be a better/fast way to do this using the java library.
            % Useage: >> timestruct = geovar.timewindowij([2004 1 1 0 0 0], [2005 12 31 0 0 0]);
            %              >> timestruct = geovar.timewindowij(731947, 732677);
            % Result: time.time = time data
            %            time.index = time indices
            s = src.size;
            first = ones(1, length(s));
            last = s;
            stride = first;
            g.time = src.gettimedata(first, last, stride);

            if isfield(g, 'time') % are any of the fields recognized as time explictly
                if nargin > 2 % If two times are input, do window
                    starttime = datenum(varargin{1});
                    stoptime = datenum(varargin{2});
                    if isempty(starttime)
                        starttime = g.time(1);
                    end
                    if isempty(stoptime)
                        stoptime = g.time(end);
                    end

                    t_index1 = g.time >= starttime;
                    t_index2 = g.time <= stoptime;
                    d.index = find(t_index1==t_index2);
                    d.time = g.time(d.index);
                elseif nargin < 3 % If theres only one time, do nearest
                    neartime = datenum(varargin{1});
                    diff = abs(g.time-neartime);
                    ind = find(diff==min(diff));
                    d.index = ind(1);
                    d.time = g.time(d.index);
                    if length(ind) > 1;
                        warning('NCGEOVARIABLE:TIMEWINDOWIJ',...
                            ['Multiple time indices determined to be nearest to the supplied time,'...
                            'only the first index was output.']);
                    end
                end
            else
                me = MException(['NCTOOLBOX:ncgeovariable:timewindowij'], ...
                    'No grid variable returned as time.');
                me.throw;
            end
        end % end timewindowij

        function d = geosubset(obj, struct, varargin)
            % NCGEOVARIABLE.GEOSUBSET -
            % Create subset structure:
            % s.time=[now-3 now]; % Can be omitted or s.t_index=[i1 i2]; can be used for subsetting time using indices
            % s.t_stride=2; % Can be omitted
            % s.lat=[22 27];
            % s.lon=[117 123];
            % s.h_stride=[2 2]; % Can be omitted
            % s.z_index=[1 30]; % Can be omitted
            % or
            % s.z_index='top'
            % s.z_index='bottom'
            % s.v_stride=2; % Can be omitted
            %
            % Usage:
            % sub = gvar.geosubset(s); % Subset method
            %
            % Returns:
            % Structure with fields 'data' for data values, and 'grid' for structure of coordinate variables
            % (that come from grid_interop).

            %           if ~regexp(obj.dataset.attribute('Conventions'), 'UGRID')
            nums = obj.size;
            if isfield(struct, 'h_stride');
                if length(struct.h_stride) < 2
                    struct.h_stride(1, 2) = struct.h_stride(1,1); % Resoltution to Issue 24
                end
            else
                struct.h_stride = [1 1];
            end

            if isfield(struct, 'v_stride');
            else
                struct.v_stride = 1;
            end

            if isfield(struct, 't_stride');
            else
                struct.t_stride = 1;
            end

            [indstart_r indend_r indstart_c indend_c] = obj.geoij(struct);

            if ~isempty(indstart_r)
                if isfield(struct, 'time') % Deal with time (values) or t_index (indices) bounds
                    if iscell(struct.time)
                        switch length(struct.time)
                            case 1
                                t = obj.timewindowij(struct.time{1});
                                tmin_i = t.index;
                                tmax_i = t.index;
                            case 2
                                t = obj.timewindowij(struct.time{1}, struct.time{2});
                                tmin_i = min(t.index);
                                tmax_i = max(t.index);
                        end
                    else
                        if ischar(struct.time)
                            t = obj.timewindowij(struct.time);
                            tmin_i = t.index;
                            tmax_i = t.index;
                        elseif length(struct.time) == 6
                            t = obj.timewindowij(struct.time);
                            tmin_i = t.index;
                            tmax_i = t.index;
                        else
                            switch length(struct.time)
                                case 1
                                    t = obj.timewindowij(struct.time);
                                    tmin_i = t.index;
                                    tmax_i = t.index;
                                case 2
                                    t = obj.timewindowij(struct.time(1), struct.time(2));
                                    tmin_i = min(t.index);
                                    tmax_i = max(t.index);
                                otherwise % for anything else assume that it is a single datevec
                                    t = obj.timewindowij(struct.time);
                                    tmin_i = t.index;
                                    tmax_i = t.index;
                            end
                        end
                    end
                elseif isfield(struct, 't_index') % check for 1 time index not the same time index twice $$$$$$$$$$$$
                    if iscell(struct.t_index)
                        switch length(t_index)
                            case 1
                                if numel(struct.t_index{1}) > 1 % check to see if someone used str or datevec by accident
                                    me = MException(['NCTOOLBOX:ncgeovariable:geosubset'], ...
                                        'Expected min time to be an index/integer.');
                                    me.throw;
                                else
                                    tmin_i = struct.t_index{1};
                                end
                            case 2
                                if numel(struct.t_index{1}) > 1 % check to see if someone used str or datevec by accident
                                    me = MException(['NCTOOLBOX:ncgeovariable:geosubset'], ...
                                        'Expected min time to be an index/integer.');
                                    me.throw;
                                else
                                    tmin_i = struct.t_index{1};
                                end
                                if numel(struct.t_index{2}) > 1 % check to see if someone used str or datevec by accident
                                    me = MException(['NCTOOLBOX:' mfilename ':geosubset'], ...
                                        'Expected max time to be an index/integer.');
                                    me.throw;
                                else
                                    tmax_i = struct.t_index{2};
                                end
                        end
                    else
                        switch length(struct.t_index)
                            case 1
                                tmin_i = struct.t_index(1);
                                tmax_i = struct.t_index(1);
                            case 2
                                tmin_i = struct.t_index(1);
                                tmax_i = struct.t_index(2);
                        end
                    end
                else
                    tmin_i = 1;
                    tmax_i = nums(1);
                end

                order = obj.getaxesorder;

                if isfield(struct, 'z_index')
                    if length(struct.z_index)==1 || ischar(struct.z_index)
                        if ischar(struct.z_index)
                           e=obj.top_bot;
                           switch struct.z_index
                                case 'top'
                                    struct.z_index=e(1);
                                    struct.z_index=e(1);

                                case 'bottom'
                                    struct.z_index=e(2);
                                    struct.z_index=e(2);                         end
                        end
                        zmin = struct.z_index;
                        zmax = struct.z_index;

                    elseif length(struct.z_index)==2
                        zmin = struct.z_index(1);
                        zmax = struct.z_index(2);
                    end
                else
                    if isfield(order, 'z')
                        zmin = 1;
                        zmax = nums(order.z);
                    else
                        order.z = [];
                        zmin = [];
                        zmax = [];
                    end
                end

                if length(nums) < 2
                    me = MException(['NCTOOLBOX:ncgeovariable:geosubset'], ...
                        ['Expected data of ', obj.name, ' to be at least rank 2.']);
                    me.throw;
                else
                    first = ones([1, length(nums)]);
                    stride = ones([1, length(nums)]);
                    last = nums;

                    if order.lon ~= order.lat
                        % pass
                    else
                        order.lat = order.lon;
                        order.lon = order.lat+1;
                    end
                    try
                        stride(order.time) = struct.t_stride;
                        first(order.time)   = tmin_i;
                        last(order.time)   = tmax_i;
                    end
                    try
                        stride(order.z) = struct.v_stride;
                        first(order.z) = zmin;
                        last(order.z)         = zmax;
                    end
                    stride(order.lon)     = struct.h_stride(2);
                    stride(order.lat)      = struct.h_stride(1);

                    first(order.lon)     = indstart_c;
                    first(order.lat)      = indstart_r;

                    last(order.lon)     = indend_c;
                    last(order.lat)      = indend_r;

                    %                 else
                    %                     me = MException(['NCTOOLBOX:ncgeovariable:geosubset'], ...
                    %                         ['Expected shape of data of ', obj.name, ' to be less than rank 5.']);
                    %                     me.throw;
                end

                % Get the corresponding data and interop grid...
                d.data = obj.data(first, last, stride);
                d.grid = obj.grid_interop(first, last, stride, varargin);
                d.start = first;
                d.stop = last;
                d.stride = stride;
                for i = 1:length(first)
                    d.indices{i} = first(i):stride(i):last(i);
                end
                %           else
                %             ugrid = ncugrid(obj); % Starting to add place holders for ugrid subsetting functionality
                %             d = ugrid.unstructuredLatLonSubset(struct);
                %           end % end of ugrid if
            else
                me = MException(['NCTOOLBOX:ncgeovariable:geosubset'], ...
                    ['No indices returned cooresponding to geographic subset.']);
                me.throw;
            end
        end % end of geosubset
        %%
        function [indstart_r indend_r indstart_c indend_c] =...
                geoij(obj, struct)
            % GEOVARIABLE.GEOIJ - Function to return start and stop indices for rows and columns if the
            % the lat/lon grids are rank 2, but if the grids are vectors the function returns start and stop
            % indices of lon and then lat (in that order).
            %
            % Usage: >> [firstrow, lastrow, firstcol, lastcol] = geoij(geovar, subsetstruct)
            %            >> [firstlon, lastlon, firstlat, lastlat] = geoij(geovar, subsetstruct) % if lat/lon are vector
            %
            %
            %             s = obj.size;
            %             first = ones(1, length(s));
            %             last = s;
            %             stride = first;
            %             g = obj.grid_interop(first, last, stride); % this is doing the whole thing
            %           h = 0;
            flag = 0;

            warning off
            switch  length(obj.size)
                case 1
                    [g.lon, typex] = obj.getlondata(1, obj.size, 1);
                    [g.lat, typey] = obj.getlatdata(1, obj.size, 1);
                case 2
                    [g.lon, typex] = obj.getlondata([1, 1], obj.size, [1, 1]);
                    [g.lat, typey] = obj.getlatdata([1, 1], obj.size, [1, 1]);
                case 3
                    [g.lon, typex] = obj.getlondata([1, 1, 1], obj.size, [1, 1, 1]);
                    [g.lat, typey] = obj.getlatdata([1, 1, 1], obj.size, [1, 1, 1]);
                case 4
                    [g.lon, typex] = obj.getlondata([1, 1, 1, 1], obj.size, [1, 1, 1, 1]);
                    [g.lat, typey] = obj.getlatdata([1, 1, 1, 1], obj.size, [1, 1, 1, 1]);
            end
            warning on
            if strcmp(typex, 'GeoX')
                griddataset = ucar.nc2.dt.grid.GridDataset.open(obj.dataset.location);
                grid = griddataset.findGridByName(obj.name);
                grid = grid.getCoordinateSystem();

                %ig.y = grid.getYHorizAxis;
                %ig.x = grid.getXHorizAxis;
                [x, y] = meshgrid(g.lon,g.lat);
                s = size(x);
                x = reshape(x, [1 numel(x)]);
                y = reshape(y, [1 numel(y)]);
                tempXY = [x; y];
                projection = grid.getProjection();
                tempLatLon = projection.projToLatLon(tempXY);

                g.lat = reshape(tempLatLon(1,:), s);
                g.lon = reshape(tempLatLon(2,:), s);
            end

            g.lon = double(g.lon);
            g.lat = double(g.lat);

            if max(g.lon) > 360
                if min(g.lon) > 0
                    g.lon(g.lon>360) = g.lon(g.lon>360)-360;
                    g.lon(g.lon>360) = g.lon(g.lon>360)-360;
                    g.lon(g.lon>360) = g.lon(g.lon>360)-360;
                    g.lon(g.lon>180) = g.lon(g.lon>180)-360;
                else
                    error('NCGEOVARIABLE:GEOIJ',...
                        'Longitude contains values that follow both -180/180 and 0/360+ conventions; can not subset.');
                end
            elseif min(g.lon) < 0
                % Do nothing, we assume that input is in -180/180 conventions. And this means that the data is too.
            elseif max(g.lon) > 180
                g.lon = g.lon - 360;
            end

            %Unpack geosubset_structure
            if isfield(struct, 'lat');
                switch length(struct.lat)
                    case 1
                        flag = 1;
                    case 2
                        north_max = struct.lat(2);
                        north_min = struct.lat(1);
                end

            else
                north_max = max(max(g.lat));
                north_min = min(min(g.lat));
            end

            if isfield(struct,'lon');
                switch length(struct.lon)
                    case 1
                        flag = 1;
                    case 2
                        east_max = struct.lon(2);
                        east_min = struct.lon(1);
                end

            else
                east_max = max(max(g.lon));
                east_min = max(min(g.lon));
            end

            switch flag
                case 0
                    if ~isvector(g.lat)
                        [indlat_l1] = ((g.lat <= north_max)); %2d
                        [indlat_l2] = ((g.lat >= north_min)); %2d
                        [indlat_r indlat_c] = find((indlat_l1&indlat_l2)); % 1d each
                        indlon_l1 = zeros(size(indlat_l1));
                        indlon_l2 = zeros(size(indlat_l1));
                        for i = 1:length(indlat_r)
                            if g.lon(indlat_r(i), indlat_c(i)) <= east_max;
                                indlon_l1(indlat_r(i), indlat_c(i)) = true;
                            end
                            if g.lon(indlat_r(i), indlat_c(i)) >= east_min;
                                indlon_l2(indlat_r(i), indlat_c(i)) = true;
                            end
                        end
                        [ind_r, ind_c] = find((indlon_l1&indlon_l2));
                        indstart_c = min(ind_c); % Out
                        indend_c = max(ind_c);
                        indstart_r = min(ind_r);
                        indend_r = max(ind_r);
                    else
                        indlat1 = (g.lat <= north_max);
                        indlat2 = (g.lat >= north_min);
                        indlat = find(indlat1&indlat2);
                        indlon1 = (g.lon <= east_max);
                        indlon2 = (g.lon >= east_min);
                        indlon = find(indlon1&indlon2);
                        % Out
                        attrs = obj.dataset.attributes;
                        key = value4key(attrs, 'CF:featureType');

                        if strcmp(key, 'timeSeries')
                            a = [];
                            for i = 1:length(indlon)
                                if ~isempty(find(indlon(i)==indlat, 1))
                                    a(length(a)+1) = indlon(i);
                                end
                            end
                            indstart_c = min(a);
                            indend_c = max(a);
                            indstart_r = min(a);
                            indend_r = max(a);
                        else
                            indstart_c = min(indlon); % testing temp, i switched lon and lat here
                            indend_c = max(indlon);
                            indstart_r = min(indlat);
                            indend_r = max(indlat);
                            %
                        end
                    end
                case 1
                    if isvector(g.lat)
                        indstart_c = near(g.lon,struct.lon);
                        indstart_r = near(g.lat,struct.lat);
                        indend_c = indstart_c;
                        indend_r = indstart_r;
                    else
                        % This way does't require the addtional step, or the additional function and is verified
                        % to return the same results as using nearxy with first output and ind2ij.
                        % Ind2ij doesn't seem to add any value to the toolbox, everything is already here to
                        % get rows and columns. The third nearxy output can remain private if it confusing to
                        % users (no mention of it in help docs).
                        [~, ~, indexes] = near_lonlat(g.lon, g.lat, struct.lon, struct.lat);
                        indstart_c = indexes(2);
                        indend_c = indexes(2);
                        indstart_r = indexes(1);
                        indend_r = indexes(1);
                    end
            end
        end % end geoij

        function order = getaxesorder(obj)
            %  NCGEOVARIABLE.getaxesorder  gets the order of each axis of the ncgeovariable

            %             permute_nums = fliplr(obj.axes_info);
            ainfo = obj.axes_info;
            siz = obj.size;
            for i = 1:length(ainfo(:, 1))
                order.(ainfo{i,1}) = find(siz==ainfo{i, 2}(1));
            end

        end

        function sref = subsref(obj,s)
	    % NCGEOVARIABLE.subsref applies various .name and () methods to the ncgeovariable
            switch s(1).type
                % Use the built-in subsref for dot notation
                case '.'
                    switch s(1).subs
                        case 'gettimedata'
                            switch length(s)
                                case 1
                                    sref = obj;
                                case 2
                                    nums = obj.size;
                                    [first last stride] = indexing(s(2).subs, double(nums));
                                    sref = obj.gettimedata(first, last, stride);
                            end
                        case 'getlondata'
                            switch length(s)
                                case 1
                                    sref = obj;
                                case 2
                                    nums = obj.size;
                                    [first last stride] = indexing(s(2).subs, double(nums));
                                    sref = obj.getlondata(first, last, stride);
                            end
                        case 'getlatdata'
                            switch length(s)
                                case 1
                                    sref = obj;
                                case 2
                                    nums = obj.size;
                                    [first last stride] = indexing(s(2).subs, double(nums));
                                    sref = obj.getlatdata(first, last, stride);
                            end
                        case 'grid_interop'
                            switch length(s)
                                case 1
                                    sref = obj;
                                case 2
                                    nums = obj.size;
                                    [first last stride] = indexing(s(2).subs, double(nums));
                                    sref = obj.grid_interop(first, last, stride);
                            end
                        case 'data'
                            nums = size(obj);
                            if ~isempty(nums)
                                switch length(s)
                                    case 1
                                        sref = obj;
                                    case 2
                                        [first last stride] = indexing(s(2).subs, double(nums));
                                        sref = obj.data(first, last, stride);
                                end

                            else
                                sref = obj.data;
                                warning(['NCTOOLBOX:ncgeovariable:subsref'], ...
                                    ['Variable "' name '" has no netcdf dimension associated with it. Errors may result from non CF compliant files.'])
                            end
                        case 'grid'
                            nums = size(obj);
                            if ~isempty(nums)
                                switch length(s)
                                    case 1
                                        sref = obj;
                                    case 2
                                        [first last stride] = indexing(s(2).subs, double(nums));
                                        sref = obj.grid(first, last, stride);
                                end

                            else
                                warning(['NCTOOLBOX:ncgeovariable:subsref'], ...
                                    ['Variable "' name '" has no netcdf dimension associated with it. Errors may result from non CF compliant files.'])
                            end
                        otherwise
                            sref = builtin('subsref',obj,s);
                    end
                case '()'
                    [first last stride] = indexing(s(1).subs, double(size(obj)));
                    sref = obj.data(first, last, stride);
                case '{}'
                    warning(['NCTOOLBOX:ncgeovariable:subsref'], ...
                        'Not a supported subscripted reference, "{}" are not permitted to call variable object methods');
            end
        end
    end % methods end

    %     methods (Access = protected)
    %
    %     end % protected methods end
end % class end

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [vertices, faces] = readSTL(filename)

[fs, vs, cout] = stlreadASCII(filename);

[v,I,J] = unique(vs,'rows');

vertices = vs(sort(I),:);
faces = zeros(size(fs));
for i = 1:size(fs,2)
    vf = vs(fs(:,i), :);
    [tf, lc] = ismember(vf,vertices, 'rows');
    faces(:,i) = lc;
end

faces = unique(faces,'rows');
dif1 = faces(:,1)-faces(:,2); dif2=faces(:,2)-faces(:,3); dif3=faces(:,3)-faces(:,1);
indDif1 = find(dif1 == 0);
indDif2 = find(dif2 == 0);
indDif3 = find(dif3 == 0);
indDif = union(indDif1, indDif2, indDif3);

ufIDX = setdiff([1:size(faces,1)], indDif);
faces = faces(ufIDX,:);

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus