text:
function fmri_FFA_localizerGLMcontrasts

% LO Localizer scans processing code
%
% USER IMPORTANT!!!: Before running this script, make sure that all loloc
% scans have the correct annotation fields and scan groups (are grouped
% correctly). Otherwise, fmri_FFA_checkDataTypesForGLM will fail.
%
% 1. Set GLM parameters for all scans
% 2. Check LO localizer experiments for:
%    a. Datatype and scan numbers
%    b. Record parfiles
%    c. Make sure scan groupings are correct
% 3. Run GLM for LO localizer experiments
% 4. Run contrasts for the LO localizer experiments
%
% Calls the functions:
% * fmri_FFA_setGLMparams (subfunction)
% * fmri_FFA_checkDataTypesForGLM
% * fmri_FFA_computeContrastsLocalizer
%
%
% DY 06/05/2008
% Modified 06/18/2008: there are calls to er_setParams and saveSessions all
% over the place, which results in the params for the MotionComp_RefScan1
% being changed for the worse. I'm not sure how to efficiently prevent this
% from happening, so instead, I let it happen, and just fix the ANNOTATION
% and PARFILE fields at the end, since accuracy here is important.
%
% Code taken from fmri_FFA_motionCompFix

clear all

% Set directory and subject list
% Alternatively: subs = {'sess1','sess2','sess3',...};
if ispc
    fmriDir = 'W:\projects\Kids\fmri\localizer\';
    codeDir = 'W:\projects\Kids\fmri\code';
else
    fmriDir = '/biac2/kgs/projects/Kids/fmri/localizer/';
    codeDir = '/biac2/kgs/projects/Kids/fmri/code/';
end
cd(fmriDir); s=dir('*0*'); subs={s.name}; addpath(codeDir);

% Start a log text file to document successes and failures in preprocessing
dateAndTime=datestr(now); dateAndTime(12)='_'; dateAndTime(15)='h';dateAndTime(18)='m';
logDir=fullfile(fmriDir,'logs');
if ~isdir(logDir), mkdir(logDir), end;
logFile = fullfile(fmriDir,'logs',['LO_Log_' dateAndTime '.txt']);
fid=fopen(logFile,'w');
[junk,logFileName]=fileparts(logFile)
scriptName = 'fmri_FFA_localizerGLMcontrasts.m';
startTime = clock;

% 1. Set GLM parameters
params=fmri_FFA_setGLMparams(fid,logFileName,scriptName); % Subfunction within this script

for ii=1:length(subs)

    % Get rid of all mrVista globals in workspace. Put this at the top,
    % just in case we break from the rest of the loop.
    clear dataTYPES mrSESSION vANATOMYPATH HOMEDIR hI newDt newScan

    thisDir = fullfile(fmriDir,subs{ii}); cd(thisDir);
    fprintf(fid,'\n ------------------------------------------ \n');
    fprintf(fid,'Processing %s\n\n',thisDir);
    fprintf('Processing %s\n\n',thisDir);

    %% Check if the subject already has been processed by script (and thus
    %% should not be re-entered to this script). If they do, skip to the
    %% next subject by using CONTINUE.
    load mrSESSION.mat
    theglm=find(strcmp({dataTYPES.name},'GLMs'));
    if(~isempty(theglm)) && (isfield(dataTYPES(theglm).eventAnalysisParams,'logFile'))
        fprintf(fid,'Error: Subject already processed (has logFile field), skipping... \n\n');
        fprintf('Error: Subject already processed (has logFile field), skipping... \n\n');
        continue
    end

    %% 2a/b/c. Find datatype number for MotionComp_RefScan1, find scan
    %% numbers for lo localizer scans, record parfiles for these scans and
    %% make sure scan groupings are correct. If everything checks out,
    %% proceed with GLM. If not, SKIP THIS SUBJECT.
    [go,s1,s2]=fmri_FFA_checkDataTypesForGLM(fid,dataTYPES);
    %clear dataTYPES mrSESSION vANATOMYPATH

    % If there was a problem in step 3, skip GLM and proceed to next
    % subject
    if go==0
        continue

    %% 3. Otherwise, run GLM (to debug: fid=1)
    elseif go==1
        try
            % Initialize a hidden view of the 'inplane' data (like invoking
            % mrVista without a GUI, dataTYPES/etc in workspace)
            hI = initHiddenInplane(s1.mc, s1.scan);
            params.parfiles={s1.parfile, s2.parfile};
            er_setParams(hI, params); % write params as eventAnalysisParams fields
            tic; [hI, lo_scan] = applyGlm(hI,s1.mc,s1.scangroup,params);
            time=toc;
            fprintf(fid,'GLM COMPLETED -- dt %d, scan %d. Time: %d m %2.2f s \n\n',...
                hI.curDataType,lo_scan,floor(time/60),mod(time,60));
            fprintf('GLM COMPLETED -- dt %d, scan %d. Time: %d m %2.2f s \n\n',...
                hI.curDataType,lo_scan,floor(time/60),mod(time,60));
            glmOK=1;
        catch
            fprintf(fid,'FAILURE TO SUCCESSFULLY RUN GLM \n');
            fprintf('FAILURE TO SUCCESSFULLY RUN GLM \n');
            theerror=lasterror;
            fprintf(fid,'%s \n\n',theerror.message);
            fprintf('%s \n\n',theerror.message);
            glmOK=0;
        end

        %% 4. Compute contrast maps
        if glmOK

            try
                fmri_FFA_computeContrastsLocalizer(hI, lo_scan)
                fprintf(fid,'Computed contrasts using fmri_FFA_computeContrastsLocalizer \n');
                fprintf('Computed contrasts using fmri_FFA_computeContrastsLocalizer \n');
            catch
                fprintf(fid,'FAILURE TO SUCCESSFULLY COMPUTE CONTRASTS \n');
                fprintf('FAILURE TO SUCCESSFULLY COMPUTE CONTRASTS \n');
                theerror=lasterror;
                fprintf(fid,'%s \n\n',theerror.message);
                fprintf('%s \n\n',theerror.message);

            end
        end
    end

    % Fix the annotation and parfile fields for the MotionComp_RefScan
    % dataTYPE, since these are undesirably changed (collateral damage).
    dataTYPES(s1.mc).eventAnalysisParams(s1.scan).annotation='loloc_run1';
    dataTYPES(s2.mc).eventAnalysisParams(s2.scan).annotation='loloc_run2';

    par1=dataTYPES(s1.mc).scanParams(s1.scan).parfile;
    par2=dataTYPES(s2.mc).scanParams(s2.scan).parfile;

    dataTYPES(s1.mc).eventAnalysisParams(s1.scan).parfiles=par1;
    dataTYPES(s2.mc).eventAnalysisParams(s2.scan).parfiles=par2;
    save(fullfile(thisDir,'mrSESSION.mat'),'mrSESSION','dataTYPES', '-append');
    clear par1 par2

end

totalTime=etime(clock,startTime);

fprintf(fid,'\n ------------------------------------------ \n');
fprintf(fid,'Total running time for script: %f minutes \n',totalTime/60);
fprintf('Total running time for script: %f minutes \n',totalTime/60);

fclose(fid); % close out the log file

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function params=fmri_FFA_setGLMparams(fid,logFileName,scriptName)

% Set GLM parameters: We set the PARAMS struct fields to the values on the
% wiki, which are further explained there. Also we print the values to the
% log. http://vpnl.stanford.edu/internal/wiki/index.php/ROIanalysis#GLMs

fprintf(fid,'\n ------------------------------------------ \n');
fprintf(fid,'\n GLM PARAMETERS \n');
fprintf(fid,' ------------------------------------------ \n');

params=er_defaultParams % Initialize param struct with defaults
p.al='params.alpha=0.0500'; eval(p.al); fprintf(fid,'%s \n',p.al);
p.at='params.ampType=''betas'''; eval(p.at); fprintf(fid,'%s \n',p.at);
p.an='params.annotation=''LO localizer Runs 1&2'''; eval(p.an); fprintf(fid,'%s \n',p.an);
p.ap='params.assignParfiles=0'; eval(p.ap); fprintf(fid,'%s \n',p.ap);
p.bp='params.bslPeriod=[-6 -5 -4 -3 -2 -1 0]'; eval(p.bp); fprintf(fid,'%s \n',p.bp);
p.dt='params.detrend=1'; eval(p.dt); fprintf(fid,'%s \n',p.dt);
p.df='params.detrendFrames=20;'; eval(p.df); fprintf(fid,'%s \n',p.df);
p.ea='params.eventAnalysis=1'; eval(p.ea); fprintf(fid,'%s \n',p.ea);
p.eb='params.eventsPerBlock=6'; eval(p.eb); fprintf(fid,'%s \n',p.eb);
p.fp='params.framePeriod=2'; eval(p.fp); fprintf(fid,'%s \n',p.fp);
p.gt='params.glmHRF=3'; eval(p.gt); fprintf(fid,'%s \n',p.gt);
p.gp='params.glmHRF_params=[6 16 1 1 6 0 24]'; eval(p.gp); fprintf(fid,'%s \n',p.gp);
p.gw='params.glmWhiten=0'; eval(p.gw); fprintf(fid,'%s \n',p.gw); % temporal autocorrelation
p.ic='params.inhomoCorrect=1'; eval(p.ic); fprintf(fid,'%s \n',p.ic);
p.lp='params.lowPassFilter=0'; eval(p.lp); fprintf(fid,'%s \n',p.lp);
p.nb='params.normBsl=1'; eval(p.nb); fprintf(fid,'%s \n',p.nb);
p.od='params.onsetDelta=0'; eval(p.od); fprintf(fid,'%s \n',p.od);
p.pf='params.parfiles='''''; eval(p.pf); fprintf(fid,'%s \n',p.pf); % fill in after checkDataTypes...
p.pp='params.peakPeriod=[4 5 6 7 8 9 10 11 12 13 14]'; eval(p.pp); fprintf(fid,'%s \n',p.pp);
p.sh='params.setHRFParams=0'; eval(p.sh); fprintf(fid,'%s \n',p.sh);
p.sc='params.snrConds=1'; eval(p.sc); fprintf(fid,'%s \n',p.sc);
p.tn='params.temporalNormalization=0'; eval(p.tn); fprintf(fid,'%s \n',p.tn);
p.tw='params.timeWindow=[-6:1:24]'; eval(p.tw); fprintf(fid,'%s \n',p.tw);
p.log=['params.logFile=''' logFileName '''']; eval(p.log); fprintf(fid,'%s \n',p.log)
p.script=['params.scriptName=''' scriptName '''']; eval(p.script); fprintf(fid,'%s \n',p.script)

% params.framePeriod = dataTYPES(dt).scanParams(scan).framePeriod;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%-------------------------------------------------------
% University of Zaragoza
% Centro Politecnico Superior
% Robotics, Perception and Real Time Group
% Author:  Javier Civera -- jcivera@unizar.es
% Date   :  2007-05-09
%-------------------------------------------------------
% Returns a camera data structure containing the calibration
%-------------------------------------------------------

function cam = initialize_cam()

% Unibrain camera, calibrated on 04/09/07
d =     0.0112;
nRows = 240;
nCols = 320;
Cx =    1.7945 / d;
Cy =    1.4433 / d;
k1=     6.333e-2;
k2=     1.390e-2;
f =     2.1735;

% % Visual Compass camera
% d =     0.0112;
% nRows = 240;
% nCols = 320;
% Cx =    1.6888/d;
% Cy =    1.4393/d;
% k1=     6.239e-2;
% k2=     1.359e-2;
% f =     2.1660;

cam.k1 =    k1;
cam.k2 =    k2;
cam.nRows = nRows;
cam.nCols = nCols;
cam.Cx =    Cx;
cam.Cy =    Cy;
cam.f =     f;
cam.dx =    d;
cam.dy =    d;
cam.model = 'two_distortion_parameters';

cam.K =     sparse( [ f/d   0     Cx;
                0  f/d    Cy;
                0    0     1] );

repo: ekfmonoslam


text:
function D=envidataread(datafile,info)
%ENVIDATAREAD Reads data of ENVI image file according to ENVI image file header information.
%D=ENVIDATAREAD(DATAFILE,INFO) provides image data (D).
%D will be in whatever number format (double, int, etc.) as in the ENVI
%file.

%Original version by Ian Howat, Ohio State Universtiy, ihowat@gmail.com
%Thanks to Yushin Ahn and Ray Jung
%Heavily modified by Felix Totir.

if nargin < 2
    error('You must specify the information about data (output of envihdrread)');
end

%-improved- assert(info.header_offset==0,'Non-zero header offset not supported');

%% Set binary format parameters
switch info.byte_order
    case {0}
        machine = 'ieee-le';
    case {1}
        machine = 'ieee-be';
    otherwise
        machine = 'n';
end

iscx=false; %if it is complex
switch info.data_type
    case {1}
        format = 'uint8';
    case {2}
        format= 'int16';
    case{3}
        format= 'int32';
    case {4}
        format= 'single';
    case {5}
        format= 'double';
    case {6}
        iscx=true;
        format= 'single';
    case {9}
        iscx=true;
        format= 'double';
    case {12}
        format= 'uint16';
    case {13}
        format= 'uint32';
    case {14}
        format= 'int64';
    case {15}
        format= 'uint64';
    otherwise
        error(['File type number: ',num2str(dtype),' not supported']);
end

% *** BSQ Format ***
% [Band 1]
% R1: C1, C2, C3, ...
% R2: C1, C2, C3, ...
%  ...
% RN: C1, C2, C3, ...
%
% [Band 2]
%  ...
% [Band N]

% *** BIL Format ***
% [Row 1]
% B1: C1, C2, C3, ...
% B2: C1, C2, C3, ...
%
%  ...
% [Row N]

% *** BIP Format ***
% Row 1
% C1: B1 B2 B3, ...
% C2: B1 B2 B3, ...
% ...
% Row N

n = info.lines*info.samples*info.bands;

fid=fopen(datafile,'r');
fread(fid,info.header_offset,'uint8',0,machine); %we skip the header
if ~iscx
    D = fread(fid,n,format,0,machine); %alternatively, multibandreader (matlab function) might be used
else
    D = fread(fid,2*n,format,0,machine); %alternatively, multibandreader (matlab function) might be used
    D = complex(D(1:2:end),D(2:2:end));
end
fclose(fid);

switch lower(info.interleave)
    case {'bsq'}
        D = reshape(D,[info.samples,info.lines,info.bands]);
        D = permute(D,[2,1,3]);
    case {'bil'}
        D = reshape(D,[info.samples,info.bands,info.lines]);
        D = permute(D,[3,1,2]);
    case {'bip'}
        D = reshape(D,[info.bands,info.samples,info.lines]);
        D=permute(D,[3,2,1]);
end

author: Sable

repo: mcbench-benchmarks


text:
%
%  Convert a column vector to the vector in the domain D
%   (the inverse of the function RepresentationVector)
%
%   Syntax:  v = RepresentationInverse(z, B)
%
%   Input:  z -- (vector, sparse) column vector in sparse format
%                    representing the vector v
%           B -- (matrix/string/cell) Basis info of the vector space:
%                (i) If B is an mxn matrix, the vector space consists of
%                    mxn matrices spanned by nonzero entries of B
%                (ii) If B is a character string, the vector space consists
%                    of polynomials spanned by the support of B
%                (iii) If B is a 1xm cell, the vector space is a
%                   a Cartesian product of m vector spaces, each component
%                   is either a matrix space or a polynomial space
%                 Specs of B:
%                 * length(B) = number of spaces in the Cartesian product
%                 * each cell entry B{j} contains one of the following
%                   (i) B{j} is an mxn matrix: the j-th component of u
%                      is an mxn matrix spanned by nonzero entries of B
%                   (ii) B{j} is a character string: the j-th component of
%                      u is a polynomial spanned by the support of B
%                   (iii)B{j} is a 1x2 cell: the j-th component of u is
%                      a polynomial with
%                             B{j}{1} = cell of variable names such as
%                                          {'x','y','z'}
%                             B{j}{2} = tuple degree such as [3 2 5]
%                                 for degrees in each variable
%                         implying the complete monomial basis of the
%                         variables under the tuple degree for u{j}
%                   (iv) a 1x3 cell with variable names, tuple degree,
%                         term indices for a polynomial component u{j}
%  Output:  v -- (matrix/string/cell) a (general) vector that is either
%                (i) an mxn matrix, or
%                (ii) a character string representing a polynomial, or
%                (iii) a cell array, each v{i} is either (i) or (ii) above
%
% Example:
%
% >> full(z).'   % ans =
%     0     2     3     1     0     4     5     0     0    -7     0     6
%
% >> u = RepresentationInverse(z, {ones(3,2),'1+x+x*y+x*y^2+x^3+x^5'});
% >> u{:}  % ans =
%      0     1
%      2     0
%      3     4
% ans =
% 5 - 7*x^5 + 6*x*y^2
%

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function divdif_test08 ( )

%*****************************************************************************80
%
%% DIVDIF_TEST08 tests R8POLY_SHIFT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 October 2006
%
%  Author:
%
%    John Burkardt
%
  n = 3;

  scale = 2.0;
  shift = +3.0;
  poly_cof(1:3) = [ +6.0, -1.0, 2.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'DIVDIF_TEST08\n' );
  fprintf ( 1, '  R8POLY_SHIFT shifts polynomial coefficients.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Polynomial coefficients for argument X\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %3d  %14f\n', i, poly_cof(i) );
  end

  poly_cof = r8poly_shift ( scale, shift, n, poly_cof );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  SCALE = %f\n', scale );
  fprintf ( 1, '  SHIFT = %f\n', shift );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Polynomial coefficients for argument\n' );
  fprintf ( 1, '    Z = SCALE * X + SHIFT\n' );
  fprintf ( 1, '\n' );

  for i = 1 : n
    fprintf ( 1, '  %3d  %14f\n', i, poly_cof(i) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function stat = memprofile(varargin)

% MEMPROFILE this is a dummy placeholder

if nargout
  stat.mem = [0 0 0 0];
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [varargout]=minesweeper(varargin)
%initialization
clc
clear all
close all
% playing a background music
% % [bk,fs]=wavread('bk2.wav',1000);
% % wavplay(bk,fs,'async');
%declaration and initializing variables
h=struct;
h.row=10;
h.column=10;
h.iposx=100;h.iposy=100;
h.boxdim=[35 35];
h.gap=[0.5 0.5];
h.mine=imread(strcat('mine','.png'));
% initializing gui
h.main=figure('name','Suleman''s minesweeper','NumberTitle','off','Position',[h.iposx,h.iposy,(h.row+4)*h.boxdim(1),(h.column+4)*h.boxdim(2)]);
h.newgame=uicontrol(h.main,'style','pushbutton','position',[0,h.boxdim(2)*(h.column+4)-35,100,35],'string','New Game','callback','minesweeper','tooltipstring','Newgame');
h.close=uicontrol(h.main,'style','pushbutton','position',[h.boxdim(1)*(h.row+4)-100,0,100,35],'string','Close','callback','close gcbf');
h.popup=uicontrol(h.main,'style','popup','position',[h.boxdim(1)*(h.row+4)-150,h.boxdim(2)*(h.column+4)-20,100,20],...
    'tooltipstring','Go for Expert','string','Difficulty|Beginner|Intermediate|Expert','interruptible','off');
for x=1:h.row+2
    for y=1:h.column+2
        h.box(y,x)=uicontrol(h.main,'style','pushbutton','FontWeight','bold','foregroundcolor','b','fontsize',12,...
            'position',[h.boxdim(1)*(x),h.boxdim(2)*(h.column+3-y),h.boxdim(1)-h.gap(1),h.boxdim(2)-h.gap(2)]);
        if x==1 || y==1 || x==h.row+2 || y==h.column+2
            set(h.box(y,x),'visible','off');
        end
    end
end
% initializing background work like bomb placement
h.difficulty=0.65;
bomb=int8(h.difficulty*rand(h.row+2,h.column+2));
bomb(:,1)=zeros(h.row+2,1);bomb(:,h.column+2)=zeros(h.row+2,1);bomb(1,:)=zeros(1,h.column+2);bomb(h.row+2,:)=zeros(1,h.column+2);
h.game=bomb*9;
for x=2:h.row+1
    for y=2:h.column+1
        n=0;
        if bomb(x,y)==0
            if bomb(x-1,y-1)==1
                n=n+1;
            end
            if bomb(x-1,y)==1
                n=n+1;
            end
            if bomb(x-1,y+1)==1
                n=n+1;
            end
            if bomb(x,y-1)==1
                n=n+1;
            end
            if bomb(x,y+1)==1
                n=n+1;
            end
            if bomb(x+1,y-1)==1
                n=n+1;
            end
            if bomb(x+1,y)==1
                n=n+1;
            end
            if bomb(x+1,y+1)==1
                n=n+1;
            end
        h.game(x,y)=n;
        end
    end
end
h.game=flipud(h.game');
% disp(h.game);
% defining callbacks
set(h.popup,'callback',{@difficultylevel,h});
for ii=1:h.row+2
    for jj=1:h.column+2
        set(h.box(ii,jj),'callback',{'button',h},'buttondownfcn',{'mark',h});
    end
end
end


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

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym atan (@var{x})
%% Symbolic atan function.
%%
%% Example:
%% @example
%% @group
%% syms x
%% y = atan (x)
%%   @result{} y = (sym) atan(x)
%% @end group
%% @end example
%%
%% Note: this file is autogenerated: if you want to edit it, you might
%% want to make changes to 'generate_functions.py' instead.
%%
%% @end defmethod

function y = atan(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('atan', x);
end

%!error atan (sym(1), 2)
%!assert (isequaln (atan (sym(nan)), sym(nan)))

%!shared x, d
%! d = 1;
%! x = sym('1');

%!test
%! f1 = atan(x);
%! f2 = atan(d);
%! assert( abs(double(f1) - f2) < 1e-15 )

%!test
%! D = [d d; d d];
%! A = [x x; x x];
%! f1 = atan(A);
%! f2 = atan(D);
%! assert( all(all( abs(double(f1) - f2) < 1e-15 )))

%!test
%! % round trip
%! y = sym('y');
%! A = atan (d);
%! f = atan (y);
%! h = function_handle (f);
%! B = h (d);
%! assert (A, B, -eps)

repo: octsympy


text:
function [cost, grad] = getCostGrad(problem, x, storedb, key)
% Computes the cost function and the gradient at x in one call if possible.
%
% function [cost, grad] = getCostGrad(problem, x)
% function [cost, grad] = getCostGrad(problem, x, storedb)
% function [cost, grad] = getCostGrad(problem, x, storedb, key)
%
% Returns the value at x of the cost function described in the problem
% structure, as well as the gradient at x.
%
% storedb is a StoreDB object, key is the StoreDB key to point x.
%
% See also: canGetCost canGetGradient getCost getGradient

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:
%
%   April 3, 2015 (NB):
%       Works with the new StoreDB class system.
%
%   Aug. 2, 2018 (NB):
%       The value of the cost function is now always cached.
%
%   Sep. 6, 2018 (NB):
%       The gradient is now also cached.

    % Allow omission of the key, and even of storedb.
    if ~exist('key', 'var')
        if ~exist('storedb', 'var')
            storedb = StoreDB();
        end
        key = storedb.getNewKey();
    end

    % Contrary to most similar functions, here, we get the store by
    % default. This is for the caching functionality described below.
    store = storedb.getWithShared(key);
    store_is_stale = false;

    % Check if the cost or gradient are readily available from the store.
    force_grad_caching = true;
    if isfield(store, 'cost__')
        cost = store.cost__;
        if force_grad_caching && isfield(store, 'grad__')
            grad = store.grad__;
            return;
        else
            grad = getGradient(problem, x, storedb, key); % caches grad
            return;
        end
    end
    % If we get here, the cost was not previously cached, but maybe the
    % gradient was?
    if force_grad_caching && isfield(store, 'grad__')
        grad = store.grad__;
        cost = getCost(problem, x, storedb, key); % this call caches cost
        return;
    end

    % Neither the cost nor the gradient were available: let's compute both.

    if isfield(problem, 'costgrad')
    %% Compute the cost/grad pair using costgrad.

        % Check whether this function wants to deal with storedb or not.
        switch nargin(problem.costgrad)
            case 1
                [cost, grad] = problem.costgrad(x);
            case 2
                [cost, grad, store] = problem.costgrad(x, store);
            case 3
                % Pass along the whole storedb (by reference), with key.
                [cost, grad] = problem.costgrad(x, storedb, key);
                store_is_stale = true;
            otherwise
                up = MException('manopt:getCostGrad:badcostgrad', ...
                    'costgrad should accept 1, 2 or 3 inputs.');
                throw(up);
        end

    else
    %% Revert to calling getCost and getGradient separately

        % The two following calls will already cache cost and grad, then
        % the caches will be overwritten at the end of this function, with
        % the same values (it is not a problem).
        cost = getCost(problem, x, storedb, key);
        grad = getGradient(problem, x, storedb, key);
        store_is_stale = true;

    end

    if store_is_stale
        store = storedb.getWithShared(key);
    end

    % Cache here.
    store.cost__ = cost;
    if force_grad_caching
        store.grad__ = grad;
    end

    storedb.setWithShared(store, key);

end

author: NicolasBoumal

repo: manopt


text:
function prob = nmix2row(prob)
%NMIX2ROW  Convert Mixed Nonlinear Constraints to Nonlinear Row bounds
%   prob = nmix2row(prob)
%
%   Use -1 for <=, 0 for =, and 1 for >= in vector e

%   Copyright (C) 2012 Jonathan Currie (I2C2)

%Assign common vars
rhs = prob.nlrhs;
e = prob.nle;

%Transpose as Required
if(size(rhs,2) > 1)
    rhs = rhs';
end
if(size(e,2) > 1)
    e = e';
end
if(xor(isempty(rhs),isempty(e)))
    error('You must supply both nlrhs and nle!');
end
if(length(rhs) ~= length(e))
    error('nlrhs and nle are not the same length!');
end
if(any(e < -1) || any(e > 1))
    error('The vector e must only contain values -1, 0 or 1');
end

%Defaults
ncon = length(e);
prob.cl = -Inf(ncon,1);
prob.cu = Inf(ncon,1);

%Indices
ieq = e == 0;
icl = (e ==  1 | ieq);
icu = (e == -1 | ieq);

%Fill In
prob.cl(icl) = rhs(icl);
prob.cu(icu) = rhs(icu);

%Remove unused fields
prob.nlrhs = []; prob.nle = [];

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function test_pull2111

% WALLTIME 00:10:00
% MEM 1gb
% DEPENDENCY ft_read_data, ft_read_header

filename = dccnpath('/home/common/matlab/fieldtrip/data/test/pull2110/data_test_pull2111.mat');

hdr = ft_read_header(filename, 'headerformat', 'neuroomega_mat', 'chantype', 'ecog');
data = ft_read_data(filename, 'header', hdr, 'dataformat', 'neuroomega_mat', 'headerformat', 'neuroomega_mat');


text:
function [G] = Gal2G(Gal)
% Convert acceleration from galileos to acceration from Earths gravity.
% Chad A. Greene 2012
G = Gal*0.0010197162;


text:
function mix = mygmminit(x,centres,idx,covartype)
%GMMINIT Initialises Gaussian mixture model from data
%
[ndata, xdim] = size(x);
k = rows(centres);
mix = gmm(xdim,k,covartype);
I = eye(k);
post = I(idx,:);
mix.centres = centres;

% Arbitrary width used if variance collapses to zero: make it 'large' so
% that centre is responsible for a reasonable number of points.
GMM_WIDTH = 1.0;

% Set priors depending on number of points in each cluster
cluster_sizes = max(sum(post, 1), 1);  % Make sure that no prior is zero
mix.priors = cluster_sizes/sum(cluster_sizes); % Normalise priors

switch mix.covar_type
case 'spherical'
   if mix.ncentres > 1
      % Determine widths as distance to nearest centre
      % (or a constant if this is zero)
      for i = 1:mix.ncentres
        mix.covars(i) = mean(dist2(x(idx==i,:),mix.centres(i,:)));
      end
   else
      % Just use variance of all data points averaged over all
      % dimensions
      mix.covars = mean(diag(cov(x)));
   end
  case 'diag'
    for j = 1:mix.ncentres
      % Pick out data points belonging to this centre
      c = x(find(post(:, j)),:);
      diffs = c - (ones(size(c, 1), 1) * mix.centres(j, :));
      mix.covars(j, :) = sum((diffs.*diffs), 1)/size(c, 1);
      % Replace small entries by GMM_WIDTH value
      mix.covars(j, :) = mix.covars(j, :) + GMM_WIDTH.*(mix.covars(j, :)<eps);
    end
  case 'full'
    for j = 1:mix.ncentres
      % Pick out data points belonging to this centre
      c = x(find(post(:, j)),:);
      diffs = c - (ones(size(c, 1), 1) * mix.centres(j, :));
      mix.covars(:,:,j) = (diffs'*diffs)/(size(c, 1));
      % Add GMM_WIDTH*Identity to rank-deficient covariance matrices
      if rank(mix.covars(:,:,j)) < mix.nin
	mix.covars(:,:,j) = mix.covars(:,:,j) + GMM_WIDTH.*eye(mix.nin);
      end
    end
  case 'ppca'
    for j = 1:mix.ncentres
      % Pick out data points belonging to this centre
      c = x(find(post(:,j)),:);
      diffs = c - (ones(size(c, 1), 1) * mix.centres(j, :));
      [tempcovars, tempU, templambda] = ...
	ppca((diffs'*diffs)/size(c, 1), mix.ppca_dim);
      if length(templambda) ~= mix.ppca_dim
	error('Unable to extract enough components');
      else
        mix.covars(j) = tempcovars;
        mix.U(:, :, j) = tempU;
        mix.lambda(j, :) = templambda;
      end
    end
  otherwise
    error(['Unknown covariance type ', mix.covar_type]);
end

author: kristinbranson

repo: JAABA


text:
function [f,handles] = TS_Plot2d(featureData,TimeSeries,featureLabels,annotateParams,...
                                            showDistr,cfnParams)
% TS_Plot2d   Plot a dataset in a two-dimensional space.
%
% e.g., The space of two chosen features, or two principal components.
%
%---INPUTS:
% featureData, an Nx2 vector of where to plot each of the N data objects in the
%           two-dimensional space
%
% TimeSeries, table of time-series metadata
%
% featureLabels, cell of labels for each feature
%
% annotateParams, a structure containing all the information about how to annotate
%           data points. Fields can include (cf. BF_AnnotatePoints):
%               - n, the number of data points to annotate
%               - userInput, 0: randomly selected datapoints, 1: user clicks to annotate datapoints
%               - fdim, 1x2 vector with width and height of time series as fraction of plot size
%               - maxL, maximum length of annotated time series
%               - textAnnotation: 'Name', 'ID', or 'none' to annotate this data
%               - cmap, a cell of colors, with elements for each group
%               - theMarkerSize, a custom marker size
%               - theLineWidth: line width for annotated time series
%
% showDistr, if true (default), plots marginal density estimates for each variable
%                   (above and to the right of the plot), otherwise set to false.
%
% cfnParams, (optional) custom classification parameters to fit.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

% ------------------------------------------------------------------------------
%% Check Inputs:
% ------------------------------------------------------------------------------

% featureData should be a Nx2 vector of where to plot each of the N data objects
% in the two-dimensional space:
if nargin < 1
    error('You must provide two-dimensional feature vectors for the data.')
end

if nargin < 3 || isempty(featureLabels)
    featureLabels = {'',''};
end

if nargin < 4 || isempty(annotateParams)
    annotateParams = struct('n',6); % annotate for six points
end

% By default, plot kernel density estimates above and on the side of the plot:
if nargin < 5 || isempty(showDistr)
    showDistr = true;
end

if nargin < 6
    cfnParams = GiveMeDefaultClassificationParams(TimeSeries);
end
if isfield(cfnParams,'numFolds')
    cfnParams.numFolds = 0;
end

if isfield(annotateParams,'makeFigure')
    makeFigure = annotateParams.makeFigure;
else
    makeFigure = true; % default is to plot on a brand new figure('color','w')
end

%-------------------------------------------------------------------------------
% Preliminaries
%-------------------------------------------------------------------------------
groupLabels = TS_ExtractGroupLabels(TimeSeries);
classLabels = categories(groupLabels);
numClasses = length(classLabels);

%-------------------------------------------------------------------------------
%% Plot
%-------------------------------------------------------------------------------
if makeFigure % can set extras.makeFigure = 0 to plot within a given setting
    f = figure('color','w');
    box('on'); % white figure
    f.Position(3:4) = [600, 550];
else
    f = gcf;
end

% Set colors
if (numClasses==0) || (numClasses==1)
    groupColors = {[0,0,0]}; % Just use black...
else
    if isstruct(annotateParams) && isfield(annotateParams,'cmap')
        if iscell(annotateParams.cmap)
            groupColors = annotateParams.cmap; % specify the cell itself
        elseif ischar(annotateParams.cmap)
            groupColors = BF_GetColorMap(annotateParams.cmap,numClasses,1);
        end
    else
        groupColors = GiveMeColors(numClasses);
    end
end
annotateParams.groupColors = groupColors;

% ------------------------------------------------------------------------------
%% Plot distributions
% ------------------------------------------------------------------------------
if showDistr
    % Top distribution (marginal of first feature)
    subplot(4,4,1:3); hold('on'); box('on')
    maxx = 0; minn = 100;
    for i = 1:numClasses
        fr = BF_plot_ks(featureData(groupLabels==classLabels{i},1),groupColors{i},0);
        maxx = max([maxx,fr]); minn = min([minn,fr]);
    end
    axTop = gca;
    set(axTop,'XTickLabel',[]);
    set(axTop,'YTickLabel',[]);
    set(axTop,'ylim',[minn,maxx]);

    % Side distribution (marginal of second feature)
    subplot(4,4,[8,12,16]); hold on; box('on')
    maxx = 0; minn = 100;
    for i = 1:numClasses
        fr = BF_plot_ks(featureData(groupLabels==classLabels{i},2),groupColors{i},1);
        maxx = max([maxx,fr]); minn = min([minn,fr]);
    end
    axSide = gca;
    set(axSide,'XTickLabel',[]);
    set(axSide,'YTickLabel',[]);
    set(axSide,'xlim',[minn,maxx]);
end

% ------------------------------------------------------------------------------
%% Set up a 2D plot
% ------------------------------------------------------------------------------
if showDistr
    subplot(4,4,[5:7,9:11,13:15]); box('on');
    axMain = gca;
end
hold('on');

if isfield(annotateParams,'theMarkerSize');
    theMarkerSize = annotateParams.theMarkerSize; % specify custom marker size
else
    theMarkerSize = 12; % Marker size for '.'
end

handles = cell(numClasses,1);
for i = 1:numClasses
    handles{i} = plot(featureData(groupLabels==classLabels{i},1),...
                featureData(groupLabels==classLabels{i},2),...
                '.','color',groupColors{i},'MarkerSize',theMarkerSize);
end

% Link axes
if showDistr
    linkaxes([axMain,axTop],'x');
    linkaxes([axMain,axSide],'y');
end

%-------------------------------------------------------------------------------
% Annotate points:
%-------------------------------------------------------------------------------
% Label axes first without classification rates so the user can see what they're doing when annotating
labelAxes(0);
if ismember('Data',TimeSeries.Properties.VariableNames)
    % Annotate if time-series data is provided

    % Produce xy points
    xy = featureData;

    % Go-go-go:
    BF_AnnotatePoints(xy,TimeSeries,annotateParams);
end

% ------------------------------------------------------------------------------
%% Do classification and plot a classify boundary?
% ------------------------------------------------------------------------------
didClassify = false;
if numClasses > 1 && ~isempty(fields(cfnParams))
    % Compute the in-sample classification rate:
    classRate = nan(3,1); % classRate1, classRate2, classRateboth
    try
        fprintf(1,'Estimating %u-class classification rates for each feature (and in combination)...\n',numClasses);
        classRate(1) = GiveMeCfn(featureData(:,1),groupLabels,featureData(:,1),groupLabels,cfnParams);
        classRate(2) = GiveMeCfn(featureData(:,2),groupLabels,featureData(:,2),groupLabels,cfnParams);
        classRate(3) = GiveMeCfn(featureData,groupLabels,featureData(:,1:2),groupLabels,cfnParams);
        % Record that classification was performed successfully:
        didClassify = true;
        fprintf(1,'%s in 2-dim space: %.2f%%\n',cfnParams.whatLoss,classRate(3));
    catch emsg
        fprintf(1,'\nClassification rates not computed\n(%s)\n',emsg.message);
    end

    % Also plot an SVM classification boundary:
    if numClasses < 5
        fprintf(1,'Estimating classification boundaries...\n');
        % Train the model (in-sample):
        [~,Mdl] = GiveMeCfn(featureData,groupLabels,featureData,groupLabels,cfnParams);

        % Predict scores over the 150x150 grid through space
        gridInc = 150;
        [x1Grid,x2Grid] = meshgrid(linspace(min(featureData(:,1)),max(featureData(:,1)),gridInc),...
                                   linspace(min(featureData(:,2)),max(featureData(:,2)),gridInc));
        fullGrid = [x1Grid(:),x2Grid(:)];
        predLabels = predict(Mdl,fullGrid);

        % For each class plot the contour of their region of the space:
        if numClasses==2
            integerPredLabels = arrayfun(@(x)find(classLabels==x),predLabels);
            contour(x1Grid,x2Grid,reshape(integerPredLabels-1.5,size(x1Grid)),[0 0],...
                                    '--k','LineWidth',2);
        else
            for i = 1:numClasses
                isMostProbable = (predLabels==classLabels{i});
                if ~all(isMostProbable==isMostProbable(1));
                    contour(x1Grid,x2Grid,reshape(isMostProbable,size(x1Grid)),...
                            [0.5 0.5],'-','LineWidth',2,'color',groupColors{i});
                end
            end
        end
    end
end

%-------------------------------------------------------------------------------
% Set Legend
%-------------------------------------------------------------------------------
if numClasses > 1
    legendText = cell(numClasses,1);
    for i = 1:numClasses
        if ~isempty(classLabels)
            legendText{i} = sprintf('%s (%u)',classLabels{i},sum(groupLabels==classLabels{i}));
        else
            legendText{i} = sprintf('Group %u (%u)',i,sum(groupLabels==classLabels{i}));
        end
    end
    legend([handles{:}],legendText,'interpreter','none','Location','best','AutoUpdate','off');
end

%-------------------------------------------------------------------------------
% Relabel axes with classification rates and set title if classification
% performed successfully:
%-------------------------------------------------------------------------------
if didClassify
    labelAxes(1);
    title(sprintf('Combined classification rate (%s) = %.2f%%',cfnParams.whatClassifier, ...
                    round(classRate(3,1))),'interpreter','none');
end

%-------------------------------------------------------------------------------
%-------------------------------------------------------------------------------
function labelAxes(didClassify)
    % Label Axes
    if didClassify
        labelText = cell(2,1);
        for i = 1:2
            labelText{i} = sprintf('%s (acc = %.2f %%)',featureLabels{i}, ...
                                    round(classRate(i,1)));
        end
    else
        labelText = featureLabels;
    end

    xlabel(labelText{1},'interpreter','none')
    ylabel(labelText{2},'interpreter','none')
end
%-------------------------------------------------------------------------------

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function range = p18_range ( )

%*****************************************************************************80
%
%% P18_RANGE returns an interval bounding the root for problem 18.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real RANGE(2), the minimum and maximum values of
%    an interval containing the root.
%
  range(1) = 0.988;
  range(2) = 1.012;

  return
end


text:
function [M0,M1,L,X,q0] = spm_mfa(M,x,u)
% Jacobian for mean field approximations
% FORMAT [M0,M1,L,X,q0] = spm_mfa(M,x,u)
%--------------------------------------------------------------------------
% M   - model specification structure
% Required fields:
%    M.f   - dx/dt    = f(x,u,P)                     {function string or m-file}
%    M.g   - y(t)     = g(x,u,P)                     {function string or m-file}
%    M.m   - m inputs
%    M.n   - n states
%    M.l   - l ouputs
%    M.x   - (n x 1) = x(0) = expansion point
%    M.W   - (n x n) - covariance matrix of deterministic noise
% x    - cell array of vectors specifying evaluation grid
% u    - expansion point for inputs (c.f. background activity);
%
% M0   - 1st order Bilinear matrix dq/dt = M0*q + u*M1*q,  q = p(X);
% M1   - 2nd order Bilinear matrix
% L    - output matrix                 <y> = L*q;
% X    - evaluation points of state space
% q0   - stable mode M0*q0 = 0
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_mfa.m 4936 2012-09-18 19:47:55Z karl $

% default expansion point for inputs
%--------------------------------------------------------------------------
if nargin < 3
    u0 = zeros(M.m,1);
else
    u0 = u;
end

% create X - coordinates of evaluation grid
%--------------------------------------------------------------------------
n     = length(x);
for i = 1:n
    q     = 1;
    for j = 1:n
        q = kron(x{j}(:).^(i == j),q);
    end
    X(:,i) = q;
end

% f(x,0)
%--------------------------------------------------------------------------
N     = length(X);
f     = sparse(n,N);
for i = 1:N
    f(:,i) = feval(M.f,X(i,:)',u0,M.pE);
end

% Jacobian J - M0
%==========================================================================
for i = 1:n
    d(i)  = length(x{i});
    dx(i) = x{i}(2) - x{i}(1);
    I{i}  = speye(d(i),d(i));
end
J     = sparse(N,N);

% cycle over dimensions of state space
%--------------------------------------------------------------------------
for i = 1:length(x)

    % differential operators (for positive and negative flows)
    %----------------------------------------------------------------------
    j          = find(f(i,:) < 0);
    u          = sparse(j,1,f(i,j),N,1);
    j          = find(f(i,:) > 0);
    v          = sparse(j,1,f(i,j),N,1);

    dp         = spdiags(ones(d(i),1)*[1 -1],[ 0;1],d(i),d(i))/dx(i);
    dp(:,1)    = 0;
    dq         = spdiags(ones(d(i),1)*[1 -1],[-1;0],d(i),d(i))/dx(i);
    dq(:,end)  = 0;

    % Kronecker tensor products
    %----------------------------------------------------------------------
    Dp         = 1;
    Dq         = 1;
    for j = 1:(i - 1)
        Dp = kron(I{j},Dp);
        Dq = kron(I{j},Dq);
    end
    Dp    = kron(dp,Dp);
    Dq    = kron(dq,Dq);
    for j = (i + 1):n
        Dp = kron(I{j},Dp);
        Dq = kron(I{j},Dq);
    end
    DP{i} = Dp;

    % augment Jacobian
    %----------------------------------------------------------------------
    J     = J + Dp*spdiags(u,0,N,N) + Dq*spdiags(v,0,N,N);

end

% dispersion
%--------------------------------------------------------------------------
for i = 1:n
    for j = 1:n
        if M.W(i,j)
            J  = J - M.W(i,j)*DP{i}*DP{j}'/2;
        end

    end
end

% return if only M0 is required
%--------------------------------------------------------------------------
M0  = J;
if nargout == 1,  return,  end

% input induced changes to Jacobian dJ/du - M1
%==========================================================================
du    = 1e-6;
M1    = {};
for i = 1:M.m
    u     = u0;
    u(i)  = du;
    Mu    = spm_mfa(M,x,u);
    M1{i} = (Mu - M0)/du;
end

% return if only M0, M1 are required
%--------------------------------------------------------------------------
if nargout == 2,  return,  end

% Output matrix L - M1
%==========================================================================

% l(x,0)
%--------------------------------------------------------------------------
L     = sparse(M.l,N);
for i = 1:N
    L(:,i) = feval(M.g,X(i,:)',u,M.pE);
end

% equilibrium mode
%--------------------------------------------------------------------------
if nargout == 5

    % stable mode - stochastic iteration
    %----------------------------------------------------------------------
    q   = ones(N,1)/N;
    for i = 1:1024
        dq = M0*q;
        q  = q + dq/N;
        q  = q/sum(q);
    end
    q0  = q;

    % analytic solution (not implemented)
    %----------------------------------------------------------------------
    % q   = length(M0);
    % J   = [M0; ones(1,q)];
    % q0  = sum(inv(J'*J),2);

end

author: spm


text:
function [scores, blocks] = readscorestats(inffile, model)

% [scores, blocks] = readscorestats(inffile, model)
%
% Read a score statistics file and parse it into labels, total scores,
% unique flags, and block scores.

fid = fopen(inffile, 'rb');
num = fread(fid, 1, 'int32');
scores = zeros(num, 1);
blocks = zeros(num, model.numblocks);
for i = 1:num
  tmp = fread(fid, model.numblocks+1, 'double');
  scores(i) = tmp(1);
  blocks(i,:) = tmp(2:end);
end
fclose(fid);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function varargout = exportfig(varargin)
%EXPORTFIG  Export a figure.
%   EXPORTFIG(H, FILENAME) writes the figure H to FILENAME.  H is
%   a figure handle and FILENAME is a string that specifies the
%   name of the output file.
%
%   EXPORTFIG(H, FILENAME, OPTIONS) writes the figure H to FILENAME
%   with options initially specified by the structure OPTIONS. The
%   field names of OPTIONS must be legal parameters listed below
%   and the field values must be legal values for the corresponding
%   parameter. Default options can be set in releases prior to R12
%   by storing the OPTIONS structure in the root object's appdata
%   with the command
%      setappdata(0,'exportfigdefaults', OPTIONS)
%   and for releases after R12 by setting the preference with the
%   command
%      setpref('exportfig', 'defaults', OPTIONS)
%
%   EXPORTFIG(...,PARAM1,VAL1,PARAM2,VAL2,...) specifies
%   parameters that control various characteristics of the output
%   file. Any parameter value can be the string 'auto' which means
%   the parameter uses the default factory behavior, overriding
%   any other default for the parameter.
%
%   Format Paramter:
%     'Format'  a string
%          specifies the output format. Defaults to 'eps'. For a
%          list of export formats type 'help print'.
%     'Preview' one of the strings 'none', 'tiff'
%          specifies a preview for EPS files. Defaults to 'none'.
%
%   Size Parameters:
%     'Width'   a positive scalar
%          specifies the width in the figure's PaperUnits
%     'Height'  a positive scalar
%          specifies the height in the figure's PaperUnits
%     'Bounds' one of the strings 'tight', 'loose'
%          specifies a tight or loose bounding box. Defaults to 'tight'.
%     'Reference' an axes handle or a string
%          specifies that the width and height parameters
%          are relative to the given axes. If a string is
%          specified then it must evaluate to an axes handle.
%
%     Specifying only one dimension sets the other dimension
%     so that the exported aspect ratio is the same as the
%     figure's or reference axes' current aspect ratio.
%     If neither dimension is specified the size defaults to
%     the width and height from the figure's or reference
%     axes' size. Tight bounding boxes are only computed for
%     2-D views and in that case the computed bounds enclose all
%     text objects.
%
%   Rendering Parameters:
%     'Color'     one of the strings 'bw', 'gray', 'cmyk'
%         'bw'    specifies that lines and text are exported in
%                 black and all other objects in grayscale
%         'gray'  specifies that all objects are exported in grayscale
%         'rgb'   specifies that all objects are exported in color
%                 using the RGB color space
%         'cmyk'  specifies that all objects are exported in color
%                 using the CMYK color space
%     'Renderer'  one of 'painters', 'zbuffer', 'opengl'
%         specifies the renderer to use
%     'Resolution'   a positive scalar
%         specifies the resolution in dots-per-inch.
%     'LockAxes'  one of 0 or 1
%         specifies that all axes limits and ticks should be fixed
%         while exporting.
%
%     The default color setting is 'bw'.
%
%   Font Parameters:
%     'FontMode'     one of the strings 'scaled', 'fixed'
%     'FontSize'     a positive scalar
%          in 'scaled' mode multiplies with the font size of each
%          text object to obtain the exported font size
%          in 'fixed' mode specifies the font size of all text
%          objects in points
%     'DefaultFixedFontSize' a positive scalar
%          in 'fixed' mode specified the default font size in
%          points
%     'FontSizeMin' a positive scalar
%          specifies the minimum font size allowed after scaling
%     'FontSizeMax' a positive scalar
%          specifies the maximum font size allowed after scaling
%     'FontEncoding' one of the strings 'latin1', 'adobe'
%          specifies the character encoding of the font
%     'SeparateText' one of 0 or 1
%          specifies that the text objects are stored in separate
%          file as EPS with the base filename having '_t' appended.
%
%     If FontMode is 'scaled' but FontSize is not specified then a
%     scaling factor is computed from the ratio of the size of the
%     exported figure to the size of the actual figure.
%
%     The default 'FontMode' setting is 'scaled'.
%
%   Line Width Parameters:
%     'LineMode'     one of the strings 'scaled', 'fixed'
%     'LineWidth'    a positive scalar
%     'DefaultFixedLineWidth' a positive scalar
%     'LineWidthMin' a positive scalar
%          specifies the minimum line width allowed after scaling
%     'LineWidthMax' a positive scalar
%          specifies the maximum line width allowed after scaling
%     The semantics of 'Line' parameters are exactly the
%     same as the corresponding 'Font' parameters, except that
%     they apply to line widths instead of font sizes.
%
%   Style Map Parameter:
%     'LineStyleMap'    one of [], 'bw', or a function name or handle
%          specifies how to map line colors to styles. An empty
%          style map means styles are not changed. The style map
%          'bw' is a built-in mapping that maps lines with the same
%          color to the same style and otherwise cycles through the
%          available styles. A user-specified map is a function
%          that takes as input a cell array of line objects and
%          outputs a cell array of line style strings. The default
%          map is [].
%
%   Examples:
%     exportfig(gcf,'fig1.eps','height',3);
%       Exports the current figure to the file named 'fig1.eps' with
%       a height of 3 inches (assuming the figure's PaperUnits is
%       inches) and an aspect ratio the same as the figure's aspect
%       ratio on screen.
%
%     opts = struct('FontMode','fixed','FontSize',10,'height',3);
%     exportfig(gcf, 'fig2.eps', opts, 'height', 5);
%       Exports the current figure to 'fig2.eps' with all
%       text in 10 point fonts and with height 5 inches.
%
%   See also PREVIEWFIG, APPLYTOFIG, RESTOREFIG, PRINT.

%  Copyright 2000 Ben Hinkle
%  Email bug reports and comments to bhinkle@mathworks.com

if (nargin < 2)
  error('Too few input arguments');
end

% exportfig(H, filename, [options,] ...)
H = varargin{1};
if ~LocalIsHG(H,'figure')
  error('First argument must be a handle to a figure.');
end
filename = varargin{2};
if ~ischar(filename)
  error('Second argument must be a string.');
end
paramPairs = {varargin{3:end}};
if nargin > 2
  if isstruct(paramPairs{1})
    pcell = LocalToCell(paramPairs{1});
    paramPairs = {pcell{:}, paramPairs{2:end}};
  end
end
verstr = version;
majorver = str2num(verstr(1));
defaults = [];
if majorver > 5
  if ispref('exportfig','defaults')
    defaults = getpref('exportfig','defaults');
  end
elseif exist('getappdata')
  defaults = getappdata(0,'exportfigdefaults');
end
if ~isempty(defaults)
  dcell = LocalToCell(defaults);
  paramPairs = {dcell{:}, paramPairs{:}};
end

% Do some validity checking on param-value pairs
if (rem(length(paramPairs),2) ~= 0)
  error(['Invalid input syntax. Optional parameters and values' ...
	 ' must be in pairs.']);
end

auto.format = 'eps';
auto.preview = 'none';
auto.width = -1;
auto.height = -1;
auto.color = 'bw';
auto.defaultfontsize=10;
auto.fontsize = -1;
auto.fontmode='scaled';
auto.fontmin = 8;
auto.fontmax = 60;
auto.defaultlinewidth = 1.0;
auto.linewidth = -1;
auto.linemode=[];
auto.linemin = 0.5;
auto.linemax = 100;
auto.fontencoding = 'latin1';
auto.renderer = [];
auto.resolution = [];
auto.stylemap = [];
auto.applystyle = 0;
auto.refobj = -1;
auto.bounds = 'tight';
explicitbounds = 0;
auto.lockaxes = 1;
auto.separatetext = 0;
opts = auto;

% Process param-value pairs
args = {};
for k = 1:2:length(paramPairs)
  param = lower(paramPairs{k});
  if ~ischar(param)
    error('Optional parameter names must be strings');
  end
  value = paramPairs{k+1};

  switch (param)
   case 'format'
    opts.format = LocalCheckAuto(lower(value),auto.format);
    if strcmp(opts.format,'preview')
      error(['Format ''preview'' no longer supported. Use PREVIEWFIG' ...
	     ' instead.']);
    end
   case 'preview'
    opts.preview = LocalCheckAuto(lower(value),auto.preview);
    if ~strcmp(opts.preview,{'none','tiff'})
      error('Preview must be ''none'' or ''tiff''.');
    end
   case 'width'
    opts.width = LocalToNum(value, auto.width);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~LocalIsPositiveScalar(opts.width)
	error('Width must be a numeric scalar > 0');
      end
    end
   case 'height'
    opts.height = LocalToNum(value, auto.height);
    if ~ischar(value) | ~strcmp(value,'auto')
      if(~LocalIsPositiveScalar(opts.height))
	error('Height must be a numeric scalar > 0');
      end
    end
   case 'color'
    opts.color = LocalCheckAuto(lower(value),auto.color);
    if ~strcmp(opts.color,{'bw','gray','rgb','cmyk'})
      error('Color must be ''bw'', ''gray'',''rgb'' or ''cmyk''.');
    end
   case 'fontmode'
    opts.fontmode = LocalCheckAuto(lower(value),auto.fontmode);
    if ~strcmp(opts.fontmode,{'scaled','fixed'})
      error('FontMode must be ''scaled'' or ''fixed''.');
    end
   case 'fontsize'
    opts.fontsize = LocalToNum(value,auto.fontsize);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~LocalIsPositiveScalar(opts.fontsize)
	error('FontSize must be a numeric scalar > 0');
      end
    end
   case 'defaultfixedfontsize'
    opts.defaultfontsize = LocalToNum(value,auto.defaultfontsize);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~LocalIsPositiveScalar(opts.defaultfontsize)
	error('DefaultFixedFontSize must be a numeric scalar > 0');
      end
    end
   case 'fontsizemin'
    opts.fontmin = LocalToNum(value,auto.fontmin);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~LocalIsPositiveScalar(opts.fontmin)
	error('FontSizeMin must be a numeric scalar > 0');
      end
    end
   case 'fontsizemax'
    opts.fontmax = LocalToNum(value,auto.fontmax);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~LocalIsPositiveScalar(opts.fontmax)
	error('FontSizeMax must be a numeric scalar > 0');
      end
    end
   case 'fontencoding'
    opts.fontencoding = LocalCheckAuto(lower(value),auto.fontencoding);
    if ~strcmp(opts.fontencoding,{'latin1','adobe'})
      error('FontEncoding must be ''latin1'' or ''adobe''.');
    end
   case 'linemode'
    opts.linemode = LocalCheckAuto(lower(value),auto.linemode);
    if ~strcmp(opts.linemode,{'scaled','fixed'})
      error('LineMode must be ''scaled'' or ''fixed''.');
    end
   case 'linewidth'
    opts.linewidth = LocalToNum(value,auto.linewidth);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~LocalIsPositiveScalar(opts.linewidth)
	error('LineWidth must be a numeric scalar > 0');
      end
    end
   case 'defaultfixedlinewidth'
    opts.defaultlinewidth = LocalToNum(value,auto.defaultlinewidth);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~LocalIsPositiveScalar(opts.defaultlinewidth)
	error(['DefaultFixedLineWidth must be a numeric scalar >' ...
	       ' 0']);
      end
    end
   case 'linewidthmin'
    opts.linemin = LocalToNum(value,auto.linemin);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~LocalIsPositiveScalar(opts.linemin)
	error('LineWidthMin must be a numeric scalar > 0');
      end
    end
   case 'linewidthmax'
    opts.linemax = LocalToNum(value,auto.linemax);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~LocalIsPositiveScalar(opts.linemax)
	error('LineWidthMax must be a numeric scalar > 0');
      end
    end
   case 'linestylemap'
    opts.stylemap = LocalCheckAuto(value,auto.stylemap);
   case 'renderer'
    opts.renderer = LocalCheckAuto(lower(value),auto.renderer);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~strcmp(opts.renderer,{'painters','zbuffer','opengl'})
	error(['Renderer must be ''painters'', ''zbuffer'' or' ...
	       ' ''opengl''.']);
      end
    end
   case 'resolution'
    opts.resolution = LocalToNum(value,auto.resolution);
    if ~ischar(value) | ~strcmp(value,'auto')
      if ~(isnumeric(value) & (prod(size(value)) == 1) & (value >= 0));
	error('Resolution must be a numeric scalar >= 0');
      end
    end
   case 'applystyle' % means to apply the options and not export
    opts.applystyle = 1;
   case 'reference'
    if ischar(value)
      if strcmp(value,'auto')
	opts.refobj = auto.refobj;
      else
	opts.refobj = eval(value);
      end
    else
      opts.refobj = value;
    end
    if ~LocalIsHG(opts.refobj,'axes')
      error('Reference object must evaluate to an axes handle.');
    end
   case 'bounds'
    opts.bounds = LocalCheckAuto(lower(value),auto.bounds);
    explicitbounds = 1;
    if ~strcmp(opts.bounds,{'tight','loose'})
      error('Bounds must be ''tight'' or ''loose''.');
    end
   case 'lockaxes'
    opts.lockaxes = LocalToNum(value,auto.lockaxes);
   case 'separatetext'
    opts.separatetext = LocalToNum(value,auto.separatetext);
   otherwise
    error(['Unrecognized option ' param '.']);
  end
end

% make sure figure is up-to-date
drawnow;

allLines  = findall(H, 'type', 'line');
allText   = findall(H, 'type', 'text');
allAxes   = findall(H, 'type', 'axes');
allImages = findall(H, 'type', 'image');
allLights = findall(H, 'type', 'light');
allPatch  = findall(H, 'type', 'patch');
allSurf   = findall(H, 'type', 'surface');
allRect   = findall(H, 'type', 'rectangle');
allFont   = [allText; allAxes];
allColor  = [allLines; allText; allAxes; allLights];
allMarker = [allLines; allPatch; allSurf];
allEdge   = [allPatch; allSurf];
allCData  = [allImages; allPatch; allSurf];

old.objs = {};
old.prop = {};
old.values = {};

% Process format
if strncmp(opts.format,'eps',3) & ~strcmp(opts.preview,'none')
  args = {args{:}, ['-' opts.preview]};
end

hadError = 0;
oldwarn = warning;
try

  % lock axes limits, ticks and labels if requested
  if opts.lockaxes
    old = LocalManualAxesMode(old, allAxes, 'TickMode');
    old = LocalManualAxesMode(old, allAxes, 'TickLabelMode');
    old = LocalManualAxesMode(old, allAxes, 'LimMode');
  end

  % Process size parameters
  figurePaperUnits = get(H, 'PaperUnits');
  oldFigureUnits = get(H, 'Units');
  oldFigPos = get(H,'Position');
  set(H, 'Units', figurePaperUnits);
  figPos = get(H,'Position');
  refsize = figPos(3:4);
  if opts.refobj ~= -1
    oldUnits = get(opts.refobj, 'Units');
    set(opts.refobj, 'Units', figurePaperUnits);
    r = get(opts.refobj, 'Position');
    refsize = r(3:4);
    set(opts.refobj, 'Units', oldUnits);
  end
  aspectRatio = refsize(1)/refsize(2);
  if (opts.width == -1) & (opts.height == -1)
    opts.width = refsize(1);
    opts.height = refsize(2);
  elseif (opts.width == -1)
    opts.width = opts.height * aspectRatio;
  elseif (opts.height == -1)
    opts.height = opts.width / aspectRatio;
  end
  wscale = opts.width/refsize(1);
  hscale = opts.height/refsize(2);
  sizescale = min(wscale,hscale);
  old = LocalPushOldData(old,H,'PaperPositionMode', ...
			 get(H,'PaperPositionMode'));
  set(H, 'PaperPositionMode', 'auto');
  newPos = [figPos(1) figPos(2)+figPos(4)*(1-hscale) ...
	    wscale*figPos(3) hscale*figPos(4)];
  set(H, 'Position', newPos);
  set(H, 'Units', oldFigureUnits);

  % process line-style map
  if ~isempty(opts.stylemap) & ~isempty(allLines)
    oldlstyle = LocalGetAsCell(allLines,'LineStyle');
    old = LocalPushOldData(old, allLines, {'LineStyle'}, ...
			   oldlstyle);
    newlstyle = oldlstyle;
    if ischar(opts.stylemap) & strcmpi(opts.stylemap,'bw')
      newlstyle = LocalMapColorToStyle(allLines);
    else
      try
	newlstyle = feval(opts.stylemap,allLines);
      catch
	warning(['Skipping stylemap. ' lasterr]);
      end
    end
    set(allLines,{'LineStyle'},newlstyle);
  end

  % Process rendering parameters
  switch (opts.color)
   case {'bw', 'gray'}
    if ~strcmp(opts.color,'bw') & strncmp(opts.format,'eps',3)
      opts.format = [opts.format 'c'];
    end
    args = {args{:}, ['-d' opts.format]};

    %compute and set gray colormap
    oldcmap = get(H,'Colormap');
    newgrays = 0.30*oldcmap(:,1) + 0.59*oldcmap(:,2) + 0.11*oldcmap(:,3);
    newcmap = [newgrays newgrays newgrays];
    old = LocalPushOldData(old, H, 'Colormap', oldcmap);
    set(H, 'Colormap', newcmap);

    %compute and set ColorSpec and CData properties
    old = LocalUpdateColors(allColor, 'color', old);
    old = LocalUpdateColors(allAxes, 'xcolor', old);
    old = LocalUpdateColors(allAxes, 'ycolor', old);
    old = LocalUpdateColors(allAxes, 'zcolor', old);
    old = LocalUpdateColors(allMarker, 'MarkerEdgeColor', old);
    old = LocalUpdateColors(allMarker, 'MarkerFaceColor', old);
    old = LocalUpdateColors(allEdge, 'EdgeColor', old);
    old = LocalUpdateColors(allEdge, 'FaceColor', old);
    old = LocalUpdateColors(allCData, 'CData', old);

   case {'rgb','cmyk'}
    if strncmp(opts.format,'eps',3)
      opts.format = [opts.format 'c'];
      args = {args{:}, ['-d' opts.format]};
      if strcmp(opts.color,'cmyk')
	args = {args{:}, '-cmyk'};
      end
    else
      args = {args{:}, ['-d' opts.format]};
    end
   otherwise
    error('Invalid Color parameter');
  end
  if (~isempty(opts.renderer))
    args = {args{:}, ['-' opts.renderer]};
  end
  if (~isempty(opts.resolution)) | ~strncmp(opts.format,'eps',3)
    if isempty(opts.resolution)
      opts.resolution = 0;
    end
    args = {args{:}, ['-r' int2str(opts.resolution)]};
  end

  % Process font parameters
  if ~isempty(opts.fontmode)
    oldfonts = LocalGetAsCell(allFont,'FontSize');
    oldfontunits = LocalGetAsCell(allFont,'FontUnits');
    set(allFont,'FontUnits','points');
    switch (opts.fontmode)
     case 'fixed'
      if (opts.fontsize == -1)
	set(allFont,'FontSize',opts.defaultfontsize);
      else
	set(allFont,'FontSize',opts.fontsize);
      end
     case 'scaled'
      if (opts.fontsize == -1)
	scale = sizescale;
      else
	scale = opts.fontsize;
      end
      newfonts = LocalScale(oldfonts,scale,opts.fontmin,opts.fontmax);
      set(allFont,{'FontSize'},newfonts);
     otherwise
      error('Invalid FontMode parameter');
    end
    old = LocalPushOldData(old, allFont, {'FontSize'}, oldfonts);
    old = LocalPushOldData(old, allFont, {'FontUnits'}, oldfontunits);
  end
  if strcmp(opts.fontencoding,'adobe') & strncmp(opts.format,'eps',3)
    args = {args{:}, '-adobecset'};
  end

  % Process line parameters
  if ~isempty(opts.linemode)
    oldlines = LocalGetAsCell(allMarker,'LineWidth');
    old = LocalPushOldData(old, allMarker, {'LineWidth'}, oldlines);
    switch (opts.linemode)
     case 'fixed'
      if (opts.linewidth == -1)
	set(allMarker,'LineWidth',opts.defaultlinewidth);
      else
	set(allMarker,'LineWidth',opts.linewidth);
      end
     case 'scaled'
      if (opts.linewidth == -1)
	scale = sizescale;
      else
	scale = opts.linewidth;
      end
      newlines = LocalScale(oldlines, scale, opts.linemin, opts.linemax);
      set(allMarker,{'LineWidth'},newlines);
    end
  end

  % adjust figure bounds to surround axes
  if strcmp(opts.bounds,'tight')
    if (~strncmp(opts.format,'eps',3) & LocalHas3DPlot(allAxes)) | ...
	  (strncmp(opts.format,'eps',3) & opts.separatetext)
      if (explicitbounds == 1)
	warning(['Cannot compute ''tight'' bounds. Using ''loose''' ...
		 ' bounds.']);
      end
      opts.bounds = 'loose';
    end
  end
  warning('off');
  if ~isempty(allAxes)
    if strncmp(opts.format,'eps',3)
      if strcmp(opts.bounds,'loose')
	args = {args{:}, '-loose'};
      end
      old = LocalPushOldData(old,H,'Position', oldFigPos);
    elseif strcmp(opts.bounds,'tight')
      oldaunits = LocalGetAsCell(allAxes,'Units');
      oldapos = LocalGetAsCell(allAxes,'Position');
      oldtunits = LocalGetAsCell(allText,'units');
      oldtpos = LocalGetAsCell(allText,'Position');
      set(allAxes,'units','points');
      apos = LocalGetAsCell(allAxes,'Position');
      oldunits = get(H,'Units');
      set(H,'units','points');
      origfr = get(H,'position');
      fr = [];
      for k=1:length(allAxes)
	if ~strcmpi(get(allAxes(k),'Tag'),'legend')
	  axesR = apos{k};
	  r = LocalAxesTightBoundingBox(axesR, allAxes(k));
	  r(1:2) = r(1:2) + axesR(1:2);
	  fr = LocalUnionRect(fr,r);
	end
      end
      if isempty(fr)
	fr = [0 0 origfr(3:4)];
      end
      for k=1:length(allAxes)
	ax = allAxes(k);
	r = apos{k};
	r(1:2) = r(1:2) - fr(1:2);
	set(ax,'Position',r);
      end
      old = LocalPushOldData(old, allAxes, {'Position'}, oldapos);
      old = LocalPushOldData(old, allText, {'Position'}, oldtpos);
      old = LocalPushOldData(old, allText, {'Units'}, oldtunits);
      old = LocalPushOldData(old, allAxes, {'Units'}, oldaunits);
      old = LocalPushOldData(old, H, 'Position', oldFigPos);
      old = LocalPushOldData(old, H, 'Units', oldFigureUnits);
      r = [origfr(1) origfr(2)+origfr(4)-fr(4) fr(3:4)];
      set(H,'Position',r);
    else
      args = {args{:}, '-loose'};
      old = LocalPushOldData(old,H,'Position', oldFigPos);
    end
  end

  % Process text in a separate file if needed
  if opts.separatetext & ~opts.applystyle
    % First hide all text and export
    oldtvis = LocalGetAsCell(allText,'visible');
    set(allText,'visible','off');
    oldax = LocalGetAsCell(allAxes,'XTickLabel',1);
    olday = LocalGetAsCell(allAxes,'YTickLabel',1);
    oldaz = LocalGetAsCell(allAxes,'ZTickLabel',1);
    null = cell(length(oldax),1);
    [null{:}] = deal([]);
    set(allAxes,{'XTickLabel'},null);
    set(allAxes,{'YTickLabel'},null);
    set(allAxes,{'ZTickLabel'},null);
    print(H, filename, args{:});
    set(allText,{'Visible'},oldtvis);
    set(allAxes,{'XTickLabel'},oldax);
    set(allAxes,{'YTickLabel'},olday);
    set(allAxes,{'ZTickLabel'},oldaz);
    % Now hide all non-text and export as eps in painters
    [path, name, ext] = fileparts(filename);
    tfile = fullfile(path,[name '_t.eps']);
    tfile2 = fullfile(path,[name '_t2.eps']);
    foundRenderer = 0;
    for k=1:length(args)
      if strncmp('-d',args{k},2)
	args{k} = '-deps';
      elseif strncmp('-zbuffer',args{k},8) | ...
	    strncmp('-opengl', args{k},6)
	args{k} = '-painters';
	foundRenderer = 1;
      end
    end
    if ~foundRenderer
      args = {args{:}, '-painters'};
    end
    allNonText = [allLines; allLights; allPatch; ...
		  allImages; allSurf; allRect];
    oldvis = LocalGetAsCell(allNonText,'visible');
    oldc = LocalGetAsCell(allAxes,'color');
    oldaxg = LocalGetAsCell(allAxes,'XGrid');
    oldayg = LocalGetAsCell(allAxes,'YGrid');
    oldazg = LocalGetAsCell(allAxes,'ZGrid');
    [null{:}] = deal('off');
    set(allAxes,{'XGrid'},null);
    set(allAxes,{'YGrid'},null);
    set(allAxes,{'ZGrid'},null);
    set(allNonText,'Visible','off');
    set(allAxes,'Color','none');
    print(H, tfile2, args{:});
    set(allNonText,{'Visible'},oldvis);
    set(allAxes,{'Color'},oldc);
    set(allAxes,{'XGrid'},oldaxg);
    set(allAxes,{'YGrid'},oldayg);
    set(allAxes,{'ZGrid'},oldazg);
    %hack up the postscript file
    fid1 = fopen(tfile,'w');
    fid2 = fopen(tfile2,'r');
    line = fgetl(fid2);
    while ischar(line)
      if strncmp(line,'%%Title',7)
	fprintf(fid1,'%s\n',['%%Title: ', tfile]);
      elseif (length(line) < 3)
	fprintf(fid1,'%s\n',line);
      elseif ~strcmp(line(end-2:end),' PR') & ...
	    ~strcmp(line(end-1:end),' L')
	fprintf(fid1,'%s\n',line);
      end
      line = fgetl(fid2);
    end
    fclose(fid1);
    fclose(fid2);
    delete(tfile2);

  elseif ~opts.applystyle
    drawnow;
    print(H, filename, args{:});
  end
  warning(oldwarn);

catch
  warning(oldwarn);
  hadError = 1;
end

% Restore figure settings
if opts.applystyle
  varargout{1} = old;
else
  for n=1:length(old.objs)
    if ~iscell(old.values{n}) & iscell(old.prop{n})
      old.values{n} = {old.values{n}};
    end
    set(old.objs{n}, old.prop{n}, old.values{n});
  end
end

if hadError
  error(deblank(lasterr));
end

%
%  Local Functions
%

function outData = LocalPushOldData(inData, objs, prop, values)
outData.objs = {objs, inData.objs{:}};
outData.prop = {prop, inData.prop{:}};
outData.values = {values, inData.values{:}};

function cellArray = LocalGetAsCell(fig,prop,allowemptycell);
cellArray = get(fig,prop);
if nargin < 3
  allowemptycell = 0;
end
if ~iscell(cellArray) & (allowemptycell | ~isempty(cellArray))
  cellArray = {cellArray};
end

function newArray = LocalScale(inArray, scale, minv, maxv)
n = length(inArray);
newArray = cell(n,1);
for k=1:n
  newArray{k} = min(maxv,max(minv,scale*inArray{k}(1)));
end

function gray = LocalMapToGray1(color)
gray = color;
if ischar(color)
  switch color(1)
   case 'y'
    color = [1 1 0];
   case 'm'
    color = [1 0 1];
   case 'c'
    color = [0 1 1];
   case 'r'
    color = [1 0 0];
   case 'g'
    color = [0 1 0];
   case 'b'
    color = [0 0 1];
   case 'w'
    color = [1 1 1];
   case 'k'
    color = [0 0 0];
  end
end
if ~ischar(color)
  gray = 0.30*color(1) + 0.59*color(2) + 0.11*color(3);
end

function newArray = LocalMapToGray(inArray);
n = length(inArray);
newArray = cell(n,1);
for k=1:n
  color = inArray{k};
  if ~isempty(color)
    color = LocalMapToGray1(color);
  end
  if isempty(color) | ischar(color)
    newArray{k} = color;
  else
    newArray{k} = [color color color];
  end
end

function newArray = LocalMapColorToStyle(inArray);
inArray = LocalGetAsCell(inArray,'Color');
n = length(inArray);
newArray = cell(n,1);
styles = {'-','--',':','-.'};
uniques = [];
nstyles = length(styles);
for k=1:n
  gray = LocalMapToGray1(inArray{k});
  if isempty(gray) | ischar(gray) | gray < .05
    newArray{k} = '-';
  else
    if ~isempty(uniques) & any(gray == uniques)
      ind = find(gray==uniques);
    else
      uniques = [uniques gray];
      ind = length(uniques);
    end
    newArray{k} = styles{mod(ind-1,nstyles)+1};
  end
end

function newArray = LocalMapCData(inArray);
n = length(inArray);
newArray = cell(n,1);
for k=1:n
  color = inArray{k};
  if (ndims(color) == 3) & isa(color,'double')
    gray = 0.30*color(:,:,1) + 0.59*color(:,:,2) + 0.11*color(:,:,3);
    color(:,:,1) = gray;
    color(:,:,2) = gray;
    color(:,:,3) = gray;
  end
  newArray{k} = color;
end

function outData = LocalUpdateColors(inArray, prop, inData)
value = LocalGetAsCell(inArray,prop);
outData.objs = {inData.objs{:}, inArray};
outData.prop = {inData.prop{:}, {prop}};
outData.values = {inData.values{:}, value};
if (~isempty(value))
  if strcmp(prop,'CData')
    value = LocalMapCData(value);
  else
    value = LocalMapToGray(value);
  end
  set(inArray,{prop},value);
end

function bool = LocalIsPositiveScalar(value)
bool = isnumeric(value) & ...
       prod(size(value)) == 1 & ...
       value > 0;

function value = LocalToNum(value,auto)
if ischar(value)
  if strcmp(value,'auto')
    value = auto;
  else
    value = str2num(value);
  end
end

%convert a struct to {field1,val1,field2,val2,...}
function c = LocalToCell(s)
f = fieldnames(s);
v = struct2cell(s);
opts = cell(2,length(f));
opts(1,:) = f;
opts(2,:) = v;
c = {opts{:}};

function c = LocalIsHG(obj,hgtype)
c = 0;
if (length(obj) == 1) & ishandle(obj)
  c = strcmp(get(obj,'type'),hgtype);
end

function c = LocalHas3DPlot(a)
zticks = LocalGetAsCell(a,'ZTickLabel');
c = 0;
for k=1:length(zticks)
  if ~isempty(zticks{k})
    c = 1;
    return;
  end
end

function r = LocalUnionRect(r1,r2)
if isempty(r1)
  r = r2;
elseif isempty(r2)
  r = r1;
elseif max(r2(3:4)) > 0
  left = min(r1(1),r2(1));
  bot = min(r1(2),r2(2));
  right = max(r1(1)+r1(3),r2(1)+r2(3));
  top = max(r1(2)+r1(4),r2(2)+r2(4));
  r = [left bot right-left top-bot];
else
  r = r1;
end

function c = LocalLabelsMatchTicks(labs,ticks)
c = 0;
try
  t1 = num2str(ticks(1));
  n = length(ticks);
  tend = num2str(ticks(n));
  c = strncmp(labs(1),t1,length(labs(1))) & ...
      strncmp(labs(n),tend,length(labs(n)));
end

function r = LocalAxesTightBoundingBox(axesR, a)
r = [];
atext = findall(a,'type','text','visible','on');
if ~isempty(atext)
  set(atext,'units','points');
  res=LocalGetAsCell(atext,'extent');
  for n=1:length(atext)
    r = LocalUnionRect(r,res{n});
  end
end
if strcmp(get(a,'visible'),'on')
  r = LocalUnionRect(r,[0 0 axesR(3:4)]);
  oldunits = get(a,'fontunits');
  set(a,'fontunits','points');
  label = text(0,0,'','parent',a,...
	       'units','points',...
	       'fontsize',get(a,'fontsize'),...
	       'fontname',get(a,'fontname'),...
	       'fontweight',get(a,'fontweight'),...
	       'fontangle',get(a,'fontangle'),...
	       'visible','off');
  fs = get(a,'fontsize');

  % handle y axis tick labels
  ry = [0 -fs/2 0 axesR(4)+fs];
  ylabs = get(a,'yticklabels');
  yticks = get(a,'ytick');
  maxw = 0;
  if ~isempty(ylabs)
    for n=1:size(ylabs,1)
      set(label,'string',ylabs(n,:));
      ext = get(label,'extent');
      maxw = max(maxw,ext(3));
    end
    if ~LocalLabelsMatchTicks(ylabs,yticks) & ...
	  strcmp(get(a,'xaxislocation'),'bottom')
      ry(4) = ry(4) + 1.5*ext(4);
    end
    if strcmp(get(a,'yaxislocation'),'left')
      ry(1) = -(maxw+5);
    else
      ry(1) = axesR(3);
    end
    ry(3) = maxw+5;
    r = LocalUnionRect(r,ry);
  end

  % handle x axis tick labels
  rx = [0 0 0 fs+5];
  xlabs = get(a,'xticklabels');
  xticks = get(a,'xtick');
  if ~isempty(xlabs)
    if strcmp(get(a,'xaxislocation'),'bottom')
      rx(2) = -(fs+5);
      if ~LocalLabelsMatchTicks(xlabs,xticks);
	rx(4) = rx(4) + 2*fs;
	rx(2) = rx(2) - 2*fs;
      end
    else
      rx(2) = axesR(4);
      % exponent is still below axes
      if ~LocalLabelsMatchTicks(xlabs,xticks);
	rx(4) = rx(4) + axesR(4) + 2*fs;
	rx(2) = -2*fs;
      end
    end
    set(label,'string',xlabs(1,:));
    ext1 = get(label,'extent');
    rx(1) = -ext1(3)/2;
    set(label,'string',xlabs(size(xlabs,1),:));
    ext2 = get(label,'extent');
    rx(3) = axesR(3) + (ext2(3) + ext1(3))/2;
    r = LocalUnionRect(r,rx);
  end
  set(a,'fontunits',oldunits);
  delete(label);
end

function c = LocalManualAxesMode(old, allAxes, base)
xs = ['X' base];
ys = ['Y' base];
zs = ['Z' base];
oldXMode = LocalGetAsCell(allAxes,xs);
oldYMode = LocalGetAsCell(allAxes,ys);
oldZMode = LocalGetAsCell(allAxes,zs);
old = LocalPushOldData(old, allAxes, {xs}, oldXMode);
old = LocalPushOldData(old, allAxes, {ys}, oldYMode);
old = LocalPushOldData(old, allAxes, {zs}, oldZMode);
set(allAxes,xs,'manual');
set(allAxes,ys,'manual');
set(allAxes,zs,'manual');
c = old;

function val = LocalCheckAuto(val, auto)
if ischar(val) & strcmp(val,'auto')
  val = auto;
end