text:
classdef StimClass < matlab.mixin.Copyable

    % SNIRF-spec class properties
    properties
        name
        data
    end

    % Non-SNIRF class properties
    properties (Access = private)
        filename
        fileformat
    end

    % Properties not part of the SNIRF spec. These parameters aren't loaded or saved to files
    properties (Access = private)
        errmargin
    end

    methods

        % -------------------------------------------------------
        function obj = StimClass(varargin)
            % Set class properties not part of the SNIRF format
            obj.fileformat = 'hdf5';
            obj.errmargin = 1e-3;

            if nargin==1
                if isa(varargin{1}, 'StimClass')
                    obj.Copy(varargin{1});
                elseif ischar(varargin{1})
                    if exist(varargin{1}, 'file')==2
                        obj.filename = varargin{1};
                        obj.Load(varargin{1});
                    else
                        obj.name = varargin{1};
                        obj.data = [];
                    end
                end
            elseif nargin==2
                if ischar(varargin{1})
                    obj.filename = varargin{1};
                    obj.Load(varargin{1}, obj.fileformat, varargin{2});
                else
                    t        = varargin{1};
                    CondName = varargin{2};
                    obj.name = CondName;
                    for ii=1:length(t)
                        obj.data(end+1,:) = [t(ii), 5, 1];
                    end
                end
            elseif nargin==3
                s        = varargin{1};
                t        = varargin{2};
                CondName = varargin{3};
                obj.name = CondName;
                k = s>0 | s==-1 | s==-2;  % Include stim marks with these values
                obj.data = [t(k), 5*ones(length(t(k)),1), s(k)];
            elseif nargin==0
                obj.name = '';
                obj.data = [];
            end

        end

        % -------------------------------------------------------
        function err = LoadHdf5(obj, fileobj, location)
            err = 0;

            % Arg 1
            if ~exist('fileobj','var') || (ischar(fileobj) && ~exist(fileobj,'file'))
                fileobj = '';
            end

            % Arg 2
            if ~exist('location', 'var') || isempty(location)
                location = '/nirs/stim1';
            elseif location(1)~='/'
                location = ['/',location];
            end

            % Error checking
            if ~isempty(fileobj) && ischar(fileobj)
                obj.filename = fileobj;
            elseif isempty(fileobj)
                fileobj = obj.filename;
            end
            if isempty(fileobj)
               err = -1;
               return;
            end

            try

                % Open group
                [gid, fid] = HDF5_GroupOpen(fileobj, location);

                % Load datasets
                obj.name   = HDF5_DatasetLoad(gid, 'name');
                obj.data   = HDF5_DatasetLoad(gid, 'data', [], '2D');
                if all(obj.data(:)==0)
                    obj.data = [];
                end

                % Close group
                HDF5_GroupClose(fileobj, gid, fid);

            catch ME

                err = -1;
                return

            end

        end

        % -------------------------------------------------------
        function SaveHdf5(obj, fileobj, location)
            if isempty(obj.data)
                obj.data = 0;
            end

            % Arg 1
            if ~exist('fileobj', 'var') || isempty(fileobj)
                error('Unable to save file. No file name given.')
            end

            % Arg 2
            if ~exist('location', 'var') || isempty(location)
                location = '/nirs/stim1';
            elseif location(1)~='/'
                location = ['/',location];
            end

            if ~exist(fileobj, 'file')
                fid = H5F.create(fileobj, 'H5F_ACC_TRUNC', 'H5P_DEFAULT', 'H5P_DEFAULT');
                H5F.close(fid);
            end
            hdf5write_safe(fileobj, [location, '/name'], obj.name);

            % Since this is a writable writeable parameter AFTER it's creation, we
            % call hdf5write_safe with the 'rw' option
            hdf5write_safe(fileobj, [location, '/data'], obj.data, 'rw:2D');
        end

        % -------------------------------------------------------
        function Update(obj, fileobj, location)
            if ~exist(fileobj, 'file')
                fid = H5F.create(fileobj, 'H5F_ACC_TRUNC', 'H5P_DEFAULT', 'H5P_DEFAULT');
                H5F.close(fid);
            end
            hdf5write_safe(fileobj, [location, '/name'], obj.name);
            hdf5write_safe(fileobj, [location, '/data'], obj.data, 'w');
        end

        % -------------------------------------------------------
        function Copy(obj, obj2)
            obj.name = obj2.name;
            obj.data = obj2.data;
        end

        % -------------------------------------------------------
        function B = eq(obj, obj2)
            B = false;
            if ~strcmp(obj.name, obj2.name)
                return;
            end

            % Dimensions matter so dimensions must equal
            if ~all(size(obj.data)==size(obj2.data))
                return;
            end

            % Now check contents
            if ~all(obj.data(:)==obj2.data(:))
                return;
            end
            B = true;
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Public Set/Get methods
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % -------------------------------------------------------
        function SetName(obj, val)
            obj.name = val;
        end

        % -------------------------------------------------------
        function val = GetName(obj)
            val = obj.name;
        end

        % -------------------------------------------------------
        function SetData(obj, val)
            obj.data = val;
        end

        % -------------------------------------------------------
        function val = GetData(obj)
            val = obj.data;
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % All other public methods
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % ----------------------------------------------------------------------------------
        function b = Exists(obj, tPt)
            b = false;
            if ~exist('tPt','var')
                return;
            end
            if isempty(obj.data)
                return;
            end
            if isempty(find( abs(obj.data(:,1)-tPt) <  obj.errmargin ))
                return;
            end
            b = true;
        end

        % ----------------------------------------------------------------------------------
        function AddStims(obj, tPts, duration, val)
            if ~exist('duration','var')
                duration = 5;
            end
            if ~exist('val','var')
                val = 1;
            end
            for ii=1:length(tPts)
                if ~obj.Exists(tPts(ii))
                    obj.data(end+1,:) = [tPts(ii), duration, val];
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function EditValue(obj, tPts, val)
            if isempty(obj.data)
                return;
            end
            if ~exist('val','var')
                val = 1;
            end
            for ii=1:length(tPts)
                k = find( abs(obj.data(:,1)-tPts(ii)) < obj.errmargin );
                if isempty(k)
                    continue;
                end
                obj.data(k,3) = val;
            end
        end

        % -------------------------------------------------------
        function [ts, v] = GetStim(obj)
            ts = [];
            v = [];
            if isempty(obj.data)
                return;
            end
            ts = obj.data(:,1);
            v = obj.data(:,3);
        end

        % ----------------------------------------------------------------------------------
        function SetTpts(obj, tPts, k)
            if ~exist('k','var')
                k = 1:size(obj.data,1);
            end
            if ~exist('tPts','var')
                return;
            end
            if length(tPts)~=1 && length(tPts)~=length(k)
                return;
            end
            obj.data(k,1) = tPts;
        end

        % -------------------------------------------------------
        function tPts = GetTpts(obj, k)
            tPts = [];
            if isempty(obj.data)
                return;
            end
            if ~exist('k','var')
                k = 1:size(obj.data,1);
            end
            tPts = obj.data(k,1);
        end

        % ----------------------------------------------------------------------------------
        function SetDuration(obj, duration, tPts)
            if isempty(obj.data)
                return;
            end
            if ~exist('duration','var')
                duration = 5;
            end
            if ~exist('tPts','var')
                tPts = obj.data(:,1);
            end
            k=[];
            for ii=1:length(tPts)
                k(ii) = find( abs(obj.data(:,1)-tPts(ii)) < obj.errmargin );
                if isempty(k)
                    continue;
                end
            end
            obj.data(k,2) = duration;
        end

        % -------------------------------------------------------
        function duration = GetDuration(obj, tPts)
            duration = [];
            if isempty(obj.data)
                return;
            end
            if ~exist('tPts','var')
                tPts = obj.data(:,1);
            end
            k = [];
            for ii=1:length(tPts)
                k(ii) = find( abs(obj.data(:,1)-tPts(ii)) < obj.errmargin );
            end
            duration = obj.data(k,2);
        end

        % ----------------------------------------------------------------------------------
        function SetValues(obj, vals, tPts)
            if isempty(obj.data)
                return;
            end
            if ~exist('vals','var')
                vals = 1;
            end
            if ~exist('tPts','var')
                tPts = obj.data(:,1);
            end
            k=[];
            for ii=1:length(tPts)
                k(ii) = find( abs(obj.data(:,1)-tPts(ii)) < obj.errmargin );
                if isempty(k)
                    continue;
                end
            end
            obj.data(k,3) = vals;
        end

        % -------------------------------------------------------
        function vals = GetValues(obj, tPts)
            vals = [];
            if isempty(obj.data)
                return;
            end
            if ~exist('tPts','var')
                tPts = obj.data(:,1);
            end
            k = [];
            for ii=1:length(tPts)
                k(ii) = find( abs(obj.data(:,1)-tPts(ii)) < obj.errmargin );
            end
            vals = obj.data(k,3);
        end

        % ----------------------------------------------------------------------------------
        function DeleteStims(obj, tPts)
            if isempty(obj.data)
                return;
            end

            % Find all stims for any conditions which match the time points and
            % delete them from data.
            k = [];
            for ii=1:length(tPts)
                k = [k, find( abs(obj.data(:,1)-tPts(ii)) < obj.errmargin )];
            end
            obj.data(k,:) = [];
        end

        % ----------------------------------------------------------------------------------
        function ToggleStims(obj, tPts)
            if isempty(obj.data)
                return;
            end

            % Find all stims for any conditions which match the time points and
            % flip it's value.
            k = [];
            for ii=1:length(tPts)
                k = [k, find( abs(obj.data(:,1)-tPts(ii)) < obj.errmargin )];
            end
            obj.data(k,3) = -1*obj.data(k,3);
        end

        % ----------------------------------------------------------------------------------
        function b = IsEmpty(obj)
            b = true;
            if isempty(obj)
                return;
            end
            if isempty(obj.name)
                return;
            end
            if isempty(obj.data)
                return;
            end
            if all(obj.data(:)==0)
                return;
            end
            b = false;
        end

        % ----------------------------------------------------------------------------------
        function nbytes = MemoryRequired(obj)
            nbytes = 0;
            if isempty(obj)
                return
            end
            nbytes = sizeof(obj.name) + sizeof(obj.data) + sizeof(obj.filename) + sizeof(obj.fileformat) + sizeof(obj.supportedFomats) + 8;
        end

    end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%  INTERNAL FUNCTION: initialize the storage of filtering
%
%  ::
%
%    [Filters,K_store,iF_store,v_store]=initialize_storage(a,P,PAI,Q,p0,...
%    exo_nbr,horizon,m,nsteps,smpl,h,store_filters)
%
%  Args:
%
%     - **a** [cell]: initial conditions of the filter for each regime
%     - **P** [cell]: initial covariance of the filter for each regime
%     - **PAI** [vector]: initial probability distribution of regimes
%     - **p0** [scalar]: number of observables
%     - **exo_nbr** [scalar]: number of exogenous
%     - **horizon** [{1}|scalar]: number of anticipated steps + 1
%     - **nsteps** [scalar]: number of forecast steps
%     - **smpl** [scalar]: number of observations
%     - **store_filters** [0|1|2|3]: 0 (no storage), 1(predicted only),
%       2(predicted and updated), 3(predicted, updated and smoothed)
%
%  Returns:
%     :
%
%     - **Filters** [struct]: structure with different fields
%     - **K_store** [cell]: Place holder for Kalman gains
%     - **iF_store** [cell]: Place holder for inverses of covariance matrices
%       of forecast errors
%     - **v_store** [cell]: Place holder for forecast errors
%
%

author: jmaih


text:
function varargout = interp_gridded(transform, val, pos_to, varargin)

% INTERP_GRIDDED computes a matrix that interpolates values that were
% observed on positions in a regular 3-D grid onto positions that are
% unstructured, e.g. the vertices of a cortical sheet.
%
% Use as
%   [val]                = interp_gridded(transform, val, pos, ...) or
%   [interpmat, distmat] = interp_gridded(transform, val, pos, ...)
% where
%   transform  homogenous coordinate transformation matrix for the volume
%   val        3-D matrix with the values in the volume
%   pos        Mx3 matrix with the vertex positions onto which the data should
%              be interpolated
%
% Optional arguments are specified in key-value pairs and can be
%    projmethod   = 'nearest', 'sphere_avg', 'sphere_weighteddistance'
%    sphereradius = number
%    distmat      = NxM matrix with precomputed distances
%    inside       = indices for inside voxels (or logical array)

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

if nargin<3
  ft_error('Not enough input arguments.');
end

% get the optional arguments
projvec      = ft_getopt(varargin, 'projvec',       1);
projweight   = ft_getopt(varargin, 'projweight',    ones(size(projvec)));
projcomb     = ft_getopt(varargin, 'projcomb',      'mean'); % or max
projthresh   = ft_getopt(varargin, 'projthresh',    []);
projmethod   = ft_getopt(varargin, 'projmethod');    % required
sphereradius = ft_getopt(varargin, 'sphereradius');  % required for some projection methods
distmat      = ft_getopt(varargin, 'distmat');       % will be computed if not present
inside       = ft_getopt(varargin, 'inside');

dim    = size(val);
dimres = svd(transform(1:3,1:3)); % to reduce the number of elements in the distance matrix
npos_to   = size(pos_to,1);
npos_from = prod(dim);

if isempty(distmat)
  % compute the distance matrix
  switch projmethod
    case 'nearest'
      % determine the nearest voxel for each vertex
      sub = round(ft_warp_apply(inv(transform), pos_to, 'homogenous'));  % express
      sub(sub(:)<1) = 1;
      sub(sub(:,1)>dim(1),1) = dim(1);
      sub(sub(:,2)>dim(2),2) = dim(2);
      sub(sub(:,3)>dim(3),3) = dim(3);
      ind = sub2ind(dim, sub(:,1), sub(:,2), sub(:,3));
      distmat = sparse(1:npos_to, ind, ones(size(ind)), npos_to, npos_from);
      if ~isempty(inside)
        % only voxels inside the brain contain a meaningful functional value
        distmat = distmat(:, inside);
      end

    case {'sphere_avg', 'sphere_weighteddistance'}
      if isempty(sphereradius)
        ft_error('sphereradius should be specified');
      end

      [X, Y, Z] = ndgrid(1:dim(1), 1:dim(2), 1:dim(3));
      pos_from  = ft_warp_apply(sparse(transform), [X(:) Y(:) Z(:)]);
      % the distance only has to be computed to voxels inside the brain
      pos_from  = pos_from(inside,:);
      npos_from = size(pos_from,1);
      % compute the distance between all voxels and each surface point
      dfromsq = sum(pos_from.^2,2); % squared distance to origin
      dtosq   = sum(pos_to.^2,  2); % squared distance to origin
      maxnpnt = double(npos_to*ceil(4/3*pi*(sphereradius/max(dimres))^3)); % initial estimate of nonzero entries
      distmat = spalloc(npos_to, npos_from, maxnpnt);
      ft_progress('init', 'textbar', 'computing distance matrix');
      for j = 1:npos_to
        ft_progress(j/npos_to);
        d   = sqrt(dfromsq + dtosq(j) - 2 * pos_from * pos_to(j,:)');
        sel = find(d<sphereradius);
        distmat(j, sel) = single(d(sel)) + eps('single');
      end
      ft_progress('close');

    case 'project'
      % do nothing I believe

    otherwise
      ft_error('unsupported projection method');
  end % case projmethod
end % if isempty distmat

%------do something with the distance matrix
switch projmethod
  case 'nearest'
    projmat         = distmat;

  case 'sphere_avg'
    projmat         = distmat;
    [ind1, ind2, d] = find(projmat);
    nnz             = full(sum(spones(projmat),2));
    for k = 1:length(ind1)
      projmat(ind1(k),ind2(k)) = 1./nnz(ind1(k));
    end

  case 'sphere_weighteddistance'
    projmat         = distmat;
    [ind1, ind2, d] = find(projmat);
    projmat         = sparse(ind1, ind2, 1./d, npos_to, npos_from);
    [ind1, ind2, d] = find(projmat);
    normnz          = sqrt(full(sum(projmat.^2, 2)));
    projmat         = sparse(ind1, ind2, d./normnz(ind1), npos_to, npos_from);

  case 'project'
      % this method is Joachim's implementation that was originally in
      % ft_sourceplot, it assumes the functional data to be defined on a
      % regular 3D grid, and that the transformation to world-space is known

      % we also need the dim
      dim = ft_getopt(varargin, 'dim');
      dat = zeros(size(pos_to,1),1);

      % convert projvec in mm to a factor, assume mean distance of 70mm
      projvec = (70-projvec)/70;
      for iproj = 1:length(projvec),
        sub = round(ft_warp_apply(inv(transform), pos_to*projvec(iproj), 'homogenous'));  % express
        sub(sub(:)<1) = 1;
        sub(sub(:,1)>dim(1),1) = dim(1);
        sub(sub(:,2)>dim(2),2) = dim(2);
        sub(sub(:,3)>dim(3),3) = dim(3);
        ind = sub2ind(dim, sub(:,1), sub(:,2), sub(:,3));
        if strcmp(projcomb,'mean')
          dat = dat + projweight(iproj) * val(ind);
        elseif strcmp(projcomb,'max')
          dat  = max([dat projweight(iproj) * val(ind)],[],2);
          tmp2 = min([dat projweight(iproj) * val(ind)],[],2);
          fi   = find(dat < max(tmp2));
          val(fi) = tmp2(fi);
        else
          ft_error('undefined method to combine projections; use cfg.projcomb= mean or max')
        end
      end
      if strcmp(projcomb,'mean'),
        dat = dat/length(projvec);
      end
      %     if ~isempty(projthresh),
      %       mm=max(abs(val(:)));
      %       maskval(abs(val) < projthresh*mm) = 0;
      %     end
      %

  otherwise
    ft_error('unsupported projection method');
end  % case projmethod

if nargout==1 && ~strcmp(projmethod, 'project')
  % return the interpolated values
  varargout{1} = projmat * val(:);
elseif nargout==1
  varargout{1} = dat;
else
  % return the interpolation and the distance matrix
  varargout{1} = projmat;
  varargout{2} = distmat;
end


text:
function [retDat] = set_y(dat,outp,ind,fInd)

%   data=set_y(DATA,Y) sets the Y matrix of data
%   data=set_y(DATA,Y,INDEXES) sets the Y matrix of data for given indexes

if nargin==2 | (isempty(ind)&isempty(fInd)),
  dat.Y=outp;
else
    if isempty(ind),
        ind = [1:size(outp,1)];
    end;
    if (nargin<4)|isempty(fInd),
        fInd = [1:size(outp,2)];
    end;
    dat.Y(ind,fInd)=outp;
end
retDat=dat;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function forwardReversibleFigures2(model,directions)
% Figures of different classes of reactions: qualitatively forward -> quantitatively reversible
%
%INPUT
%model.DrGt0

%model.S
%
% directions        subsets of qualtiatively forward  -> quantiatively reversible
%   .forwardReversible
%   .forwardReversible_bydGt0
%   .forwardReversible_bydGt0LHS
%   .forwardReversible_bydGt0Mid
%   .forwardReversible_bydGt0RHS
%   .forwardReversible_byConc_zero_fixed_DrG0
%   .forwardReversible_byConc_negative_fixed_DrG0
%   .forwardReversible_byConc_positive_fixed_DrG0
%   .forwardReversible_byConc_negative_uncertain_DrG0
%   .forwardReversible_byConc_positive_uncertain_DrG0

% thorStandard

[nMet,nRxn]=size(model.S);

% if ~isfield(model,'DrGt0Mean')
%     model.DrGt0Mean=(model.DrGt0Min+model.DrGt0Max)/2;
% end
% if ~isfield(model,'DrGtMean')
%     model.DrGtMean=(model.DrGtMin+model.DrGtMax)/2;
% end

directions=model.directions;

% close all
figureMaster=1;
figure1=1;
figure2=1;
figure345=1;
figure6=1;
figure7=1;

if figureMaster
%     %make the master plot of all 7 regions, 2, 4,6 are shaded
%     X1=1:nRxn;%nnz(directions.forwardReversible);
%     %dGrt0
%     Y0=(model.dGt0Min+model.dGt0Max)/2;
%     L0=Y0-model.dGt0Min;
%     U0=model.dGt0Max-Y0;
%     %dGrt
%     Y=(model.dGtMin+model.dGtMax)/2;
%     L=Y-model.dGtMin;
%     U=model.dGtMax-Y;
%     %find the amount of reactions with normal cumulative distribution over
%     %range of dGt0
%     P = normcdf(0,Y0,L0);

    %Y0=model.DrGt0; %model.DrGt0 = model.DrGt0 + delta_pH + delta_chi;
    Y0=(model.DrGt0Min+model.DrGt0Max)/2;%old vonB11
    L0=Y0-model.DrGt0Min;
    U0=model.DrGt0Max-Y0;

    Y=model.DrGtMean; %model.DrGtMean=(model.DrGtMax+model.DrGtMin)/2;
    L=Y-model.DrGtMin;
    U=model.DrGtMax-Y;

    P=directions.forwardProbability;
    %sort by probability that a reaction is forward (puts any NaN first)
    [tmp,xip]=sort(P,'descend');

    %     only take the indices of the problematic reactions, but be sure to
    %     take them in order of descending P
    xip2=zeros(nnz(directions.forwardReversible),1);
    p=1;
    for n=1:nRxn
         if directions.forwardReversible(xip(n))
            xip2(p)=xip(n);
            p=p+1;
        end
    end
    xip=xip2;
    X1=1:length(xip);

    %replace the NaN due to zero st dev
    %NaNpLHS=nnz(directions.forwardReversible_byConc_negative_uncertain_DrG0);
    %nNaNpRHS=nnz(directions.forwardReversible_byConc_positive_uncertain_DrG0);
    nNaNpLHS=nnz(isnan(directions.forwardReversible_byConc_negative_uncertain_DrG0));
    nNaNpRHS=nnz(isnan(directions.forwardReversible_byConc_positive_uncertain_DrG0));
    if (nNaNpLHS+nNaNpRHS)~=nnz(isnan((P(directions.forwardReversible))))
        warning('A:B','Extra category of NaN P(Delta_{r}G^{primem}<0) not taken into account');
        %nans are first in the ordering of indexes
        NaNPInd=xip(1:nNaNpLHS+nNaNpRHS);
        %sorts indices of the zero std dev met by their mean dG0t
        [tmp,xipNaNPInd]=sort(model.DrGt0Mean(NaNPInd));
        %new ordering
        xip=[NaNPInd(xipNaNPInd(1:nNaNpLHS)); xip(nNaNpLHS+nNaNpRHS+1:end); NaNPInd(xipNaNPInd(nNaNpLHS+1:nNaNpLHS+nNaNpRHS))];
    end

    figure1 = figure('PaperSize',[11 8.5],'PaperOrientation','landscape');
    % Create axes
    axes1 = axes('Parent',figure1,'Color',[0.702 0.7804 1]);
    hold on;
    %upper and lower Y
    minY=min(model.DrGtMin(directions.forwardReversible));
    maxY=max(model.DrGtMax(directions.forwardReversible));
    %baselines
    PreversibleBar_byConcLHS=ones(1,nRxn)*minY;
    PreversibleBar_byConcRHS=ones(1,nRxn)*minY;
    PreversibleBar_bydGt0=ones(1,nRxn)*minY;
    %bar for 2 & 6
    PreversibleBar_byConcLHS(directions.forwardReversible_byConc_negative_uncertain_DrG0)=maxY;
    PreversibleBar_byConcRHS(directions.forwardReversible_byConc_positive_uncertain_DrG0)=maxY;
    %bar for 4
    PreversibleBar_bydGt0(directions.forwardReversible_bydGt0Mid)=maxY;
    bar_handle4=bar(X1,PreversibleBar_bydGt0(xip),1,'BaseValue',minY,'FaceColor',[0.86 0.86 0.86],'EdgeColor','none');
    bar_handle2=bar(X1,PreversibleBar_byConcLHS(xip),1,'BaseValue',minY,'FaceColor',[0.86 0.86 0.86],'EdgeColor','none');
    bar_handle6=bar(X1,PreversibleBar_byConcRHS(xip),1,'BaseValue',minY,'FaceColor',[0.86 0.86 0.86],'EdgeColor','none');

    %dGrt errorbar
    hE=errorbar(X1,Y(xip),L(xip),U(xip),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','r');
    if 0
        % adjust error bar width
        %hE_c=get(hE, 'Children'); % deprecated since 2014b
        %errorbarXData= get(hE_c(2), 'XData');
        errorbarXData=get(hE, 'XData');
        errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0;
        errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0;
        errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0;
        errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0;
        %set(hE_c(2), 'XData', errorbarXData);
        set(hE, 'XData', errorbarXData);
    end
    %dGrt0 errorbar on top and inside dGrt
    hE2=errorbar(X1,Y0(xip),L0(xip),U0(xip),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','b');
    if 0
        % adjust error bar width
        %hE_c=get(hE2, 'Children');
        %errorbarXData= get(hE_c(2), 'XData');
        errorbarXData= get(hE2, 'XData');
        errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0;
        errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0;
        errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0;
        errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0;
        %set(hE_c(2), 'XData', errorbarXData);
        set(hE2, 'XData', errorbarXData);
    end

    %mean dGrt0
    plot(X1,Y0(xip),'.','LineStyle','none','Color',[0.3412 0.7961 0.1922]);
    %zero line
    %cumulative probability that reaction is really forward, assuming a
    %normal distribution about the mean dGt0
    [AX,H1,H2]=plotyy(X1,zeros(1,length(X1)),X1,P(xip));
    set(AX(1),'YTickMode','manual','YTick',floor(minY/100)*200:100:maxY);%,'TickDirMode','manual','TickDir','out');
    set(AX(2),'YTick',[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1]);
    set(H1,'LineStyle','none')
    set(H2,'LineStyle','-','LineWidth',2,'Color','k')
    plot(X1,zeros(1,length(X1)),'w','LineWidth',2,'LineStyle','--');
    %axis limits
    %     axis(AX(1),[0 length(X1) minY maxY])
    axis(AX(1),[0 length(X1) -500 500]);
    set(AX(1),'FontSize',16)
    set(AX(2),'FontSize',16,'YColor','k')
    axis(AX(2),[0 length(X1) 0 1])
    title('Qualitatively forward, but quantitatively reversible using group contribution estimates of \Delta_{r}G^{\primem}.','FontSize',16)
    set(get(AX(1),'Ylabel'),'String','\Delta_{r}G^{\primem}  (blue)    or     \Delta_{r}G^{\prime} (red) (kJ/mol)')
    set(get(AX(2),'Ylabel'),'String','P(\Delta_{r}G^{\primem}<0)')
    set(get(AX(1),'Ylabel'),'FontSize',16)
    set(get(AX(2),'Ylabel'),'FontSize',16)
    xlabel('Reactions, sorted by \Delta_{r}G^{\primem} or P(\Delta_{r}G^{\primem}<0)');
    saveas(figure1 ,'fwdReversible','fig');
end

%qualitatively forward reactions that are quantitatively
%reversible by concentration alone (no dGt0 error)
% fprintf('%i%s\n',nnz(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorLHS),' qualitatively forward reactions that are GC quantitatively forward by dGr0t, but reversible by concentration alone (No error in GC dGr0t).');
%if figure1 && any(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorLHS)
if ishandle(figure1) && any(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorLHS)

    %make the master plot of all 7 regions, 2, 4,6 are shaded
    X1=1:nRxn;%nnz(directions.forwardReversible);
    %dGrt0
    Y0=(model.dGt0Min+model.dGt0Max)/2;
    L0=Y0-model.dGt0Min;
    U0=model.dGt0Max-Y0;
    %dGrt
    Y=(model.dGtMin+model.dGtMax)/2;
    L=Y-model.dGtMin;
    U=model.dGtMax-Y;
    %find the amount of reactions with normal cumulative distribution over
    %range of dGt0
    P = normcdf(0,Y0,L0);
    %sort by probability that a reaction is forward (puts any NaN first)
    [tmp,xip]=sort(P,'descend');

    %     only take the indices of the problematic reactions, but be sure to
    %     take them in order of descending P
    xip2=zeros(nnz(directions.forwardReversible),1);
    p=1;
    for n=1:nRxn
        if directions.forwardReversible(xip(n))
            xip2(p)=xip(n);
            p=p+1;
        end
    end
    xip=xip2;
    X1=1:length(xip);

    %replace the NaN due to zero st dev
    nNaNpLHS=nnz(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorLHS);
    nNaNpRHS=nnz(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorRHS);
    if (nNaNpLHS+nNaNpRHS)~=nnz(isnan((P(directions.forwardReversible))))
        warning('ExtraCategory');
    end
    %nans are first in the ordering of indexes
    NaNPInd=xip(1:nNaNpLHS+nNaNpRHS);
    %sorts indices of the zero std dev met by their mean dG0t
    [tmp,xipNaNPInd]=sort(Y0(NaNPInd));
    %new ordering
    xip=[NaNPInd(xipNaNPInd(1:nNaNpLHS)); xip(nNaNpLHS+nNaNpRHS+1:end); NaNPInd(xipNaNPInd(nNaNpLHS+1:nNaNpLHS+nNaNpRHS))];
    %     xip=xip(nNaNpLHS+nNaNpRHS+1:end);

    %reactions that cannot be assigned directionality
    %cuttoff for probabilities: must be reflective about 0.5;
%     dGfGCforwardReversibleBool_bydGt0Mid=P<0.6 & P>0.4 & directions.forwardReversible;

    %     only take the indices of the problematic reactions, but be sure to
    %     take them in order of descending P
    xip2=zeros(nnz(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorLHS),1);
    p=1;
    for n=1:length(xip)
        if dGfGCforwardReversibleBool_byConc_No_dGt0ErrorLHS(xip(n))
            xip2(p)=xip(n);
            p=p+1;
        end
    end
    xip=xip2;
    X1=1:length(xip);

    figure1 = figure('PaperSize',[11 8.5],'PaperOrientation','landscape');
    % Create axes
    axes1 = axes('Parent',figure1,'Color',[0.702 0.7804 1]);
    hold on;
    %upper and lower Y
    minY=min(model.dGtMin(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorLHS));
    maxY=max(model.dGtMax(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorLHS));
    %baselines
    PreversibleBar_byConcLHS=ones(1,nRxn)*minY;
    PreversibleBar_byConcRHS=ones(1,nRxn)*minY;
    PreversibleBar_bydGt0=ones(1,nRxn)*minY;
    %bar for 2 & 6
    PreversibleBar_byConcLHS(dGfGCforwardReversibleBool_byConcLHS)=maxY;
    bar_handle2=bar(X1,PreversibleBar_byConcLHS(xip),1,'BaseValue',minY,'FaceColor',[0.86 0.86 0.86],'EdgeColor','none');

    %dGrt errorbar
    hE=errorbar(X1,Y(xip),L(xip),U(xip),'LineStyle','none','LineWidth',9,'DisplayName','forwardReversible','Color','r');
    % adjust error bar width
    hE_c=get(hE, 'Children');
    errorbarXData= get(hE_c(2), 'XData');
    errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0;
    errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0;
    set(hE_c(2), 'XData', errorbarXData);
    %dGrt0 errorbar on top and inside dGrt
    hE2=errorbar(X1,Y0(xip),L0(xip),U0(xip),'LineStyle','none','LineWidth',9,'DisplayName','forwardReversible','Color','b');
    % adjust error bar width
    hE_c=get(hE2, 'Children');
    errorbarXData= get(hE_c(2), 'XData');
    errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0;
    errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0;
    set(hE_c(2), 'XData', errorbarXData);

    %mean dGrt0
    plot(X1,Y0(xip),'.','MarkerSize',20,'LineStyle','none','Color',[0.3412 0.7961 0.1922]);
    %zero line
    %cumulative probability that reaction is really forward, assuming a
    %normal distribution about the mean dGt0
    [AX,H1,H2]=plotyy(X1,zeros(1,length(X1)),X1,P(xip));
    set(AX(1),'YTickMode','manual','YTick',floor(minY/100)*200:20:maxY);%,'TickDirMode','manual','TickDir','out');
    set(AX(2),'YTick',[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1]);
    set(H1,'LineStyle','none')
    set(H2,'LineStyle','-','LineWidth',2,'Color','k')
    plot(X1,zeros(1,length(X1)),'w','LineWidth',2,'LineStyle','--');
    %axis limits
%     axis(AX(1),[0 length(X1) minY maxY])
    axis(AX(1),[0 length(X1) -200 200]);
    set(AX(1),'FontSize',16)
    set(AX(2),'FontSize',16,'YColor','k')
    axis(AX(2),[0 length(X1) 0 1])
    title({'Qualitatively forward, but quantitatively reversible. Exact negative \Delta_{r}G^{\primem},';...
            'but with a quantitatively reversible concentration range.'},'FontSize',16)
    set(get(AX(1),'Ylabel'),'String','\Delta_{r}G^{\primem}  (blue)    or     \Delta_{r}G^{\prime} (red) (kJ/mol)')
    set(get(AX(2),'Ylabel'),'String','P(\Delta_{r}G^{\primem}<0)')
    set(get(AX(1),'Ylabel'),'FontSize',16)
    set(get(AX(2),'Ylabel'),'FontSize',16)
    xlabel('Reactions, sorted by P(\Delta_{r}G^{\primem}<0)');
end

%qualitatively reverse reactions that are quantitatively
%reversible by concentration alone (no dGt0 error)
% fprintf('%i%s\n',nnz(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorRHS),' qualitatively forward reactions that are GC quantitatively reverse by dGr0t, but reversible by concentration.(No error in GC dGr0t).');
if figure7 && any(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorRHS)
    X1=1:nRxn;
    %dGrt0
    Y0=(model.dGt0Min+model.dGt0Max)/2;
    L0=Y0-model.dGt0Min;
    U0=model.dGt0Max-Y0;
    %dGrt
    Y=(model.dGtMin+model.dGtMax)/2;
    L=Y-model.dGtMin;
    U=model.dGtMax-Y;
    %find the amount of reactions with normal cumulative distribution over
    %range of dGt0
    P = normcdf(0,Y0,L0);
    %sort by probability that a reaction is forward (puts any NaN first)
    [tmp,xip]=sort(P,'descend');

    %     only take the indices of the problematic reactions, but be sure to
    %     take them in order of descending P
    xip2=zeros(nnz(directions.forwardReversible),1);
    p=1;
    for n=1:nRxn
        if directions.forwardReversible(xip(n))
            xip2(p)=xip(n);
            p=p+1;
        end
    end
    xip=xip2;
    X1=1:length(xip);

    %replace the NaN due to zero st dev
    nNaNpLHS=nnz(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorLHS);
    nNaNpRHS=nnz(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorRHS);
    if (nNaNpLHS+nNaNpRHS)~=nnz(isnan((P(directions.forwardReversible))))
        warning('ExtraCategory');
    end
    %nans are first in the ordering of indexes
    NaNPInd=xip(1:nNaNpLHS+nNaNpRHS);
    %sorts indices of the zero std dev met by their mean dG0t
    [tmp,xipNaNPInd]=sort(Y0(NaNPInd));
    %new ordering
    xip=[NaNPInd(xipNaNPInd(1:nNaNpLHS)); xip(nNaNpLHS+nNaNpRHS+1:end); NaNPInd(xipNaNPInd(nNaNpLHS+1:nNaNpLHS+nNaNpRHS))];
    %     xip=xip(nNaNpLHS+nNaNpRHS+1:end);

    %reactions that cannot be assigned directionality
    %cuttoff for probabilities: must be reflective about 0.5;
%     dGfGCforwardReversibleBool_bydGt0Mid=P<0.6 & P>0.4 & directions.forwardReversible;

    %     only take the indices of the problematic reactions, but be sure to
    %     take them in order of descending P
    xip2=zeros(nnz(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorRHS),1);
    p=1;
    for n=1:length(xip)
        if dGfGCforwardReversibleBool_byConc_No_dGt0ErrorRHS(xip(n))
            xip2(p)=xip(n);
            p=p+1;
        end
    end
    xip=xip2;
    X1=1:length(xip);

    figure1 = figure('PaperSize',[11 8.5],'PaperOrientation','landscape');
    % Create axes
    axes1 = axes('Parent',figure1,'Color',[0.702 0.7804 1]);
    hold on;
    %upper and lower Y
    minY=min(model.dGtMin(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorRHS));
    maxY=max(model.dGtMax(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorRHS));
    %baselines
    PreversibleBar_byConcLHS=ones(1,nRxn)*minY;
    PreversibleBar_byConcRHS=ones(1,nRxn)*minY;
    PreversibleBar_bydGt0=ones(1,nRxn)*minY;
    %bar for 2 & 6
    PreversibleBar_byConcLHS(dGfGCforwardReversibleBool_byConcLHS)=maxY;
    bar_handle2=bar(X1,PreversibleBar_byConcLHS(xip),1,'BaseValue',minY,'FaceColor',[0.86 0.86 0.86],'EdgeColor','none');

    %dGrt errorbar
    hE=errorbar(X1,Y(xip),L(xip),U(xip),'LineStyle','none','LineWidth',16,'DisplayName','forwardReversible','Color','r');
    % adjust error bar width
    hE_c=get(hE, 'Children');
    errorbarXData= get(hE_c(2), 'XData');
    errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0;
    errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0;
    set(hE_c(2), 'XData', errorbarXData);
    %dGrt0 errorbar on top and inside dGrt
    hE2=errorbar(X1,Y0(xip),L0(xip),U0(xip),'LineStyle','none','LineWidth',16,'DisplayName','forwardReversible','Color','b');
    % adjust error bar width
    hE_c=get(hE2, 'Children');
    errorbarXData= get(hE_c(2), 'XData');
    errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0;
    errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0;
    set(hE_c(2), 'XData', errorbarXData);

    %mean dGrt0
    plot(X1,Y0(xip),'.','MarkerSize',24,'LineStyle','none','Color',[0.3412 0.7961 0.1922]);
    %zero line
    %cumulative probability that reaction is really forward, assuming a
    %normal distribution about the mean dGt0
    [AX,H1,H2]=plotyy(X1,zeros(1,length(X1)),X1,P(xip));
    set(AX(1),'YTickMode','manual','YTick',floor(minY/100)*200:20:maxY);%,'TickDirMode','manual','TickDir','out');
    set(AX(2),'YTick',[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1]);
    set(H1,'LineStyle','none')
    set(H2,'LineStyle','-','LineWidth',2,'Color','k')
    plot(X1,zeros(1,length(X1)),'w','LineWidth',2,'LineStyle','--');
    %axis limits
    axis(AX(1),[0 (length(X1)+1) minY maxY])
%     axis(AX(1),[0 length(X1) -500 500]);
    set(AX(1),'FontSize',16)
    set(AX(2),'FontSize',16,'YColor','k')
    axis(AX(2),[0 (length(X1)+1) 0 1])
    title({'Qualitatively forward, but quantitatively reversible. Exact positive \Delta_{r}G^{\primem},';...
            'but with a quantitatively reversible concentration range.'},'FontSize',16)
    set(get(AX(1),'Ylabel'),'String','\Delta_{r}G^{\primem}  (blue)    or     \Delta_{r}G^{\prime} (red) (kJ/mol)')
    set(get(AX(2),'Ylabel'),'String','P(\Delta_{r}G^{\primem}<0)')
    set(get(AX(1),'Ylabel'),'FontSize',16)
    set(get(AX(2),'Ylabel'),'FontSize',16)
    xlabel('Reactions, sorted by P(\Delta_{r}G^{\primem}<0)');
end

%qualitatively reverse reactions that are quantitatively
%reversible by concentration alone (with dGt0 error)
% fprintf('%i%s\n',nnz(dGfGCforwardReversibleBool_byConcLHS),' qualitatively forward reactions that are GC quantitatively reverse by dGr0t, but reversible by concentration.');
if figure2 && any(dGfGCforwardReversibleBool_byConcLHS)
    %make the master plot of all 7 regions, 2, 4,6 are shaded
    X1=1:nRxn;%nnz(directions.forwardReversible);
    %dGrt0
    Y0=(model.dGt0Min+model.dGt0Max)/2;
    L0=Y0-model.dGt0Min;
    U0=model.dGt0Max-Y0;
    %dGrt
    Y=(model.dGtMin+model.dGtMax)/2;
    L=Y-model.dGtMin;
    U=model.dGtMax-Y;
    %find the amount of reactions with normal cumulative distribution over
    %range of dGt0
    P = normcdf(0,Y0,L0);
    %sort by probability that a reaction is forward (puts any NaN first)
    [tmp,xip]=sort(P,'descend');

    %     only take the indices of the problematic reactions, but be sure to
    %     take them in order of descending P
    xip2=zeros(nnz(directions.forwardReversible),1);
    p=1;
    for n=1:nRxn
        if directions.forwardReversible(xip(n))
            xip2(p)=xip(n);
            p=p+1;
        end
    end
    xip=xip2;
    X1=1:length(xip);

    %replace the NaN due to zero st dev
    nNaNpLHS=nnz(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorLHS);
    nNaNpRHS=nnz(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorRHS);
    if (nNaNpLHS+nNaNpRHS)~=nnz(isnan((P(directions.forwardReversible))))
        warning('ExtraCategory');
    end
    %nans are first in the ordering of indexes
    NaNPInd=xip(1:nNaNpLHS+nNaNpRHS);
    %sorts indices of the zero std dev met by their mean dG0t
    [tmp,xipNaNPInd]=sort(Y0(NaNPInd));
    %new ordering
    xip=[NaNPInd(xipNaNPInd(1:nNaNpLHS)); xip(nNaNpLHS+nNaNpRHS+1:end); NaNPInd(xipNaNPInd(nNaNpLHS+1:nNaNpLHS+nNaNpRHS))];
    %     xip=xip(nNaNpLHS+nNaNpRHS+1:end);

    %reactions that cannot be assigned directionality
    %cuttoff for probabilities: must be reflective about 0.5;
%     dGfGCforwardReversibleBool_bydGt0Mid=P<0.6 & P>0.4 & directions.forwardReversible;

        %     only take the indices of the problematic reactions, but be sure to
    %     take them in order of descending P
    xip2=zeros(nnz(dGfGCforwardReversibleBool_byConcLHS),1);
    p=1;
    for n=1:length(xip)
        if dGfGCforwardReversibleBool_byConcLHS(xip(n))
            xip2(p)=xip(n);
            p=p+1;
        end
    end
    xip=xip2;
    X1=1:length(xip);

    figure1 = figure('PaperSize',[11 8.5],'PaperOrientation','landscape');
    % Create axes
    axes1 = axes('Parent',figure1,'Color',[0.702 0.7804 1]);
    hold on;
    %upper and lower Y
    minY=min(model.dGtMin(dGfGCforwardReversibleBool_byConcLHS));
    maxY=max(model.dGtMax(dGfGCforwardReversibleBool_byConcLHS));
    %baselines
    PreversibleBar_byConcLHS=ones(1,nRxn)*minY;
    PreversibleBar_byConcRHS=ones(1,nRxn)*minY;
    PreversibleBar_bydGt0=ones(1,nRxn)*minY;
    %bar for 2 & 6
    PreversibleBar_byConcLHS(dGfGCforwardReversibleBool_byConcLHS)=maxY;
    PreversibleBar_byConcRHS(dGfGCforwardReversibleBool_byConcRHS)=maxY;
    %bar for 4
    PreversibleBar_bydGt0(dGfGCforwardReversibleBool_bydGt0Mid)=maxY;
    bar_handle2=bar(X1,PreversibleBar_byConcLHS(xip),1,'BaseValue',minY,'FaceColor',[0.86 0.86 0.86],'EdgeColor','none');
%     bar_handle4=bar(X1,PreversibleBar_bydGt0(xip),1,'BaseValue',minY,'FaceColor',[0.86 0.86 0.86],'EdgeColor','none');
%     bar_handle6=bar(X1,PreversibleBar_byConcRHS(xip),1,'BaseValue',minY,'FaceColor',[0.86 0.86 0.86],'EdgeColor','none');
%
    %dGrt errorbar
    hE=errorbar(X1,Y(xip),L(xip),U(xip),'LineStyle','none','LineWidth',4.5,'DisplayName','forwardReversible','Color','r');
    % adjust error bar width
    hE_c=get(hE, 'Children');
    errorbarXData= get(hE_c(2), 'XData');
    errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0;
    errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0;
    set(hE_c(2), 'XData', errorbarXData);
    %dGrt0 errorbar on top and inside dGrt
    hE2=errorbar(X1,Y0(xip),L0(xip),U0(xip),'LineStyle','none','LineWidth',4.5,'DisplayName','forwardReversible','Color','b');
    % adjust error bar width
    hE_c=get(hE2, 'Children');
    errorbarXData= get(hE_c(2), 'XData');
    errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0;
    errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0;
    set(hE_c(2), 'XData', errorbarXData);

    %mean dGrt0
    plot(X1,Y0(xip),'.','MarkerSize',9,'LineStyle','none','Color',[0.3412 0.7961 0.1922]);
    %zero line
    %cumulative probability that reaction is really forward, assuming a
    %normal distribution about the mean dGt0
    [AX,H1,H2]=plotyy(X1,zeros(1,length(X1)),X1,P(xip));
    set(AX(1),'YTickMode','manual','YTick',floor(minY/100)*200:50:maxY);%,'TickDirMode','manual','TickDir','out');
    set(AX(2),'YTick',[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1]);
    set(H1,'LineStyle','none');
    set(H2,'LineStyle','-','LineWidth',2,'Color','k');
    plot(X1,zeros(1,length(X1)),'w','LineWidth',2,'LineStyle','--');
    %axis limits
    if length(X1)~=1
        axis(AX(1),[1 length(X1) minY maxY])
    else
        axis(AX(1),[1-0.5 length(X1)+0.5 minY maxY])
    end
%     axis(AX(1),[1 length(X1) -500 500]);
    set(AX(1),'FontSize',16)
    set(AX(2),'FontSize',16,'YColor','k')
    if length(X1)~=1
        axis(AX(2),[1 length(X1) 0 1])
    else
        axis(AX(2),[1-0.5 length(X1)+0.5 0 1])
    end
    title({'Qualitatively forward, but quantitatively reversible. Negative group contribution \Delta_{r}G^{\primem} estimate,';...
            'even with uncertainty, but with a quantitatively reversible concentration range.'},'FontSize',16)
    set(get(AX(1),'Ylabel'),'String','\Delta_{r}G^{\primem}  (blue)    or     \Delta_{r}G^{\prime} (red) (kJ/mol)')
    set(get(AX(2),'Ylabel'),'String','P(\Delta_{r}G^{\primem}<0)')
    set(get(AX(1),'Ylabel'),'FontSize',16)
    set(get(AX(2),'Ylabel'),'FontSize',16)
    xlabel('Reactions, sorted by P(\Delta_{r}G^{\primem}<0)');
end

%qualitatively reverse reactions that are quantitatively
%reversible by concentration alone (with dGt0 error)
% fprintf('%i%s\n',nnz(dGfGCforwardReversibleBool_byConcRHS),' qualitatively forward reactions that are GC quantitatively reverse by dGr0t, but reversible by concentration.');
if figure6 && any(dGfGCforwardReversibleBool_byConcRHS)
    %make the master plot of all 7 regions, 2, 4,6 are shaded
    X1=1:nRxn;%nnz(directions.forwardReversible);
    %dGrt0
    Y0=(model.dGt0Min+model.dGt0Max)/2;
    L0=Y0-model.dGt0Min;
    U0=model.dGt0Max-Y0;
    %dGrt
    Y=(model.dGtMin+model.dGtMax)/2;
    L=Y-model.dGtMin;
    U=model.dGtMax-Y;
    %find the amount of reactions with normal cumulative distribution over
    %range of dGt0
    P = normcdf(0,Y0,L0);
    %sort by probability that a reaction is forward (puts any NaN first)
    [tmp,xip]=sort(P,'descend');

    %     only take the indices of the problematic reactions, but be sure to
    %     take them in order of descending P
    xip2=zeros(nnz(directions.forwardReversible),1);
    p=1;
    for n=1:nRxn
        if directions.forwardReversible(xip(n))
            xip2(p)=xip(n);
            p=p+1;
        end
    end
    xip=xip2;
    X1=1:length(xip);

    %replace the NaN due to zero st dev
    nNaNpLHS=nnz(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorLHS);
    nNaNpRHS=nnz(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorRHS);
    if (nNaNpLHS+nNaNpRHS)~=nnz(isnan((P(directions.forwardReversible))))
        warning('ExtraCategory');
    end
    %nans are first in the ordering of indexes
    NaNPInd=xip(1:nNaNpLHS+nNaNpRHS);
    %sorts indices of the zero std dev met by their mean dG0t
    [tmp,xipNaNPInd]=sort(Y0(NaNPInd));
    %new ordering
    xip=[NaNPInd(xipNaNPInd(1:nNaNpLHS)); xip(nNaNpLHS+nNaNpRHS+1:end); NaNPInd(xipNaNPInd(nNaNpLHS+1:nNaNpLHS+nNaNpRHS))];
    %     xip=xip(nNaNpLHS+nNaNpRHS+1:end);

    %reactions that cannot be assigned directionality
    %cuttoff for probabilities: must be reflective about 0.5;
%     dGfGCforwardReversibleBool_bydGt0Mid=P<0.6 & P>0.4 & directions.forwardReversible;

        %     only take the indices of the problematic reactions, but be sure to
    %     take them in order of descending P
    xip2=zeros(nnz(dGfGCforwardReversibleBool_byConcRHS),1);
    p=1;
    for n=1:length(xip)
        if dGfGCforwardReversibleBool_byConcRHS(xip(n))
            xip2(p)=xip(n);
            p=p+1;
        end
    end
    xip=xip2;
    X1=1:length(xip);

    figure1 = figure('PaperSize',[11 8.5],'PaperOrientation','landscape');
    % Create axes
    axes1 = axes('Parent',figure1,'Color',[0.702 0.7804 1]);
    hold on;
    %upper and lower Y
    minY=min(model.dGtMin(dGfGCforwardReversibleBool_byConcRHS));
    maxY=max(model.dGtMax(dGfGCforwardReversibleBool_byConcRHS));
    %baselines
    PreversibleBar_byConcLHS=ones(1,nRxn)*minY;
    PreversibleBar_byConcRHS=ones(1,nRxn)*minY;
    PreversibleBar_bydGt0=ones(1,nRxn)*minY;
    %bar for 2 & 6
    PreversibleBar_byConcLHS(dGfGCforwardReversibleBool_byConcLHS)=maxY;
    PreversibleBar_byConcRHS(dGfGCforwardReversibleBool_byConcRHS)=maxY;
    %bar for 4
    PreversibleBar_bydGt0(dGfGCforwardReversibleBool_bydGt0Mid)=maxY;
%     bar_handle2=bar(X1,PreversibleBar_byConcLHS(xip),1,'BaseValue',minY,'FaceColor',[0.86 0.86 0.86],'EdgeColor','none');
%     bar_handle4=bar(X1,PreversibleBar_bydGt0(xip),1,'BaseValue',minY,'FaceColor',[0.86 0.86 0.86],'EdgeColor','none');
    bar_handle6=bar(X1,PreversibleBar_byConcRHS(xip),1,'BaseValue',minY,'FaceColor',[0.86 0.86 0.86],'EdgeColor','none');
%
    %dGrt errorbar
    hE=errorbar(X1,Y(xip),L(xip),U(xip),'LineStyle','none','LineWidth',10,'DisplayName','forwardReversible','Color','r');
    % adjust error bar width
    hE_c=get(hE, 'Children');
    errorbarXData= get(hE_c(2), 'XData');
    errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0;
    errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0;
    set(hE_c(2), 'XData', errorbarXData);
    %dGrt0 errorbar on top and inside dGrt
    hE2=errorbar(X1,Y0(xip),L0(xip),U0(xip),'LineStyle','none','LineWidth',10,'DisplayName','forwardReversible','Color','b');
    % adjust error bar width
    hE_c=get(hE2, 'Children');
    errorbarXData= get(hE_c(2), 'XData');
    errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0;
    errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0;
    set(hE_c(2), 'XData', errorbarXData);

    %mean dGrt0
    plot(X1,Y0(xip),'.','MarkerSize',16,'LineStyle','none','Color',[0.3412 0.7961 0.1922]);
    %zero line
    %cumulative probability that reaction is really forward, assuming a
    %normal distribution about the mean dGt0
    [AX,H1,H2]=plotyy(X1,zeros(1,length(X1)),X1,P(xip));
    set(AX(1),'YTickMode','manual','YTick',floor(minY/100)*200:20:maxY);%,'TickDirMode','manual','TickDir','out');
    set(AX(2),'YTick',[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1]);
    set(H1,'LineStyle','none');
    set(H2,'LineStyle','-','LineWidth',2,'Color','k');
    plot(X1,zeros(1,length(X1)),'w','MarkerSize',10,'LineWidth',2,'LineStyle','--');
    %axis limits
    axis(AX(1),[0.48 (length(X1)+0.48) minY maxY])
%     axis(AX(1),[0.48 (length(X1)+0.48) -500 500]);
    set(AX(1),'FontSize',16)
    set(AX(2),'FontSize',16,'YColor','k')
    axis(AX(2),[0.48 (length(X1)+0.48) 0 1])
    title({'Qualitatively forward, but quantitatively reversible. Positive group contribution \Delta_{r}G^{\primem} estimate,';...
            'even with uncertainty, but with a quantitatively reversible concentration range.'},'FontSize',16)
    set(get(AX(1),'Ylabel'),'String','\Delta_{r}G^{\primem}  (blue)    or     \Delta_{r}G^{\prime} (red) (kJ/mol)')
    set(get(AX(2),'Ylabel'),'String','P(\Delta_{r}G^{\primem}<0)')
    set(get(AX(1),'Ylabel'),'FontSize',16)
    set(get(AX(2),'Ylabel'),'FontSize',16)
    xlabel('Reactions, sorted by P(\Delta_{r}G^{\primem}<0)');
end

%qualitatively forward reactions that are quantitatively reversible by
%the range of dGt0
% fprintf('%i%s\n',nnz(dGfGCforwardReversibleBool_bydGt0),' qualitatively forward reactions that are GC quantitatively reversible by range of dGt0.');
if figure345 && any(dGfGCforwardReversibleBool_bydGt0)
    %make the master plot of all 7 regions, 2, 4,6 are shaded
    X1=1:nRxn;%nnz(directions.forwardReversible);
    %dGrt0
    Y0=(model.dGt0Min+model.dGt0Max)/2;
    L0=Y0-model.dGt0Min;
    U0=model.dGt0Max-Y0;
    %dGrt
    Y=(model.dGtMin+model.dGtMax)/2;
    L=Y-model.dGtMin;
    U=model.dGtMax-Y;
    %find the amount of reactions with normal cumulative distribution over
    %range of dGt0
    P = normcdf(0,Y0,L0);
    %sort by probability that a reaction is forward (puts any NaN first)
    [tmp,xip]=sort(P,'descend');

    %     only take the indices of the problematic reactions, but be sure to
    %     take them in order of descending P
    xip2=zeros(nnz(directions.forwardReversible),1);
    p=1;
    for n=1:nRxn
        if directions.forwardReversible(xip(n))
            xip2(p)=xip(n);
            p=p+1;
        end
    end
    xip=xip2;
    X1=1:length(xip);

    %replace the NaN due to zero st dev
    nNaNpLHS=nnz(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorLHS);
    nNaNpRHS=nnz(dGfGCforwardReversibleBool_byConc_No_dGt0ErrorRHS);
    if (nNaNpLHS+nNaNpRHS)~=nnz(isnan((P(directions.forwardReversible))))
        warning('ExtraCategory');
    end
    %nans are first in the ordering of indexes
    NaNPInd=xip(1:nNaNpLHS+nNaNpRHS);
    %sorts indices of the zero std dev met by their mean dG0t
    [tmp,xipNaNPInd]=sort(Y0(NaNPInd));
    %new ordering
    xip=[NaNPInd(xipNaNPInd(1:nNaNpLHS)); xip(nNaNpLHS+nNaNpRHS+1:end); NaNPInd(xipNaNPInd(nNaNpLHS+1:nNaNpLHS+nNaNpRHS))];
    %     xip=xip(nNaNpLHS+nNaNpRHS+1:end);

    %reactions that cannot be assigned directionality
    %cuttoff for probabilities: must be reflective about 0.5;
%     dGfGCforwardReversibleBool_bydGt0Mid=P<0.6 & P>0.4 & directions.forwardReversible;

            %     only take the indices of the problematic reactions, but be sure to
    %     take them in order of descending P
    xip2=zeros(nnz(dGfGCforwardReversibleBool_bydGt0),1);% dGfGCforwardReversibleBool_byConcRHS),1);
    p=1;
    for n=1:length(xip)
        if dGfGCforwardReversibleBool_bydGt0(xip(n))
            xip2(p)=xip(n);
            p=p+1;
        end
    end
    xip=xip2;
    X1=1:length(xip);

    figure1 = figure('PaperSize',[11 8.5],'PaperOrientation','landscape');
    % Create axes
    axes1 = axes('Parent',figure1,'Color',[0.702 0.7804 1]);
    hold on;
    %upper and lower Y
    minY=min(model.dGtMin(dGfGCforwardReversibleBool_bydGt0));
    maxY=max(model.dGtMax(dGfGCforwardReversibleBool_bydGt0));

    %baselines
    PreversibleBar_byConcLHS=ones(1,nRxn)*minY;
    PreversibleBar_byConcRHS=ones(1,nRxn)*minY;
    PreversibleBar_bydGt0=ones(1,nRxn)*minY;
    %bar for 2 & 6
    PreversibleBar_byConcLHS(dGfGCforwardReversibleBool_byConcLHS)=maxY;
    PreversibleBar_byConcRHS(dGfGCforwardReversibleBool_byConcRHS)=maxY;
    %bar for 4
    PreversibleBar_bydGt0(dGfGCforwardReversibleBool_bydGt0Mid)=maxY;
    bar_handle2=bar(X1,PreversibleBar_byConcLHS(xip),1,'BaseValue',minY,'FaceColor',[0.86 0.86 0.86],'EdgeColor','none');
    bar_handle4=bar(X1,PreversibleBar_bydGt0(xip),1,'BaseValue',minY,'FaceColor',[0.86 0.86 0.86],'EdgeColor','none');
    bar_handle6=bar(X1,PreversibleBar_byConcRHS(xip),1,'BaseValue',minY,'FaceColor',[0.86 0.86 0.86],'EdgeColor','none');

    %dGrt errorbar
    hE=errorbar(X1,Y(xip),L(xip),U(xip),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','r');
    % adjust error bar width
    %hE_c=get(hE, 'Children');
    %errorbarXData= get(hE_c(2), 'XData');
    errorbarXData= get(hE, 'XData');
    errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0;
    errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0;
    %set(hE_c(2), 'XData', errorbarXData);
    set(hE, 'XData', errorbarXData);
    %dGrt0 errorbar on top and inside dGrt
    hE2=errorbar(X1,Y0(xip),L0(xip),U0(xip),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','b');
    % adjust error bar width
    %hE_c=get(hE2, 'Children');
    %errorbarXData= get(hE_c(2), 'XData');
    errorbarXData= get(hE2, 'XData');
    errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0;
    errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0;
    errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0;
    %set(hE_c(2), 'XData', errorbarXData);
    set(hE2, 'XData', errorbarXData);

    %mean dGrt0
    plot(X1,Y0(xip),'.','MarkerSize',6,'LineStyle','none','Color',[0.3412 0.7961 0.1922]);
    %zero line
    %cumulative probability that reaction is really forward, assuming a
    %normal distribution about the mean dGt0
    [AX,H1,H2]=plotyy(X1,zeros(1,length(X1)),X1,P(xip));
    set(AX(1),'YTickMode','manual','YTick',floor(minY/100)*200:100:maxY);%,'TickDirMode','manual','TickDir','out');
    set(AX(2),'YTick',[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1]);
    set(H1,'LineStyle','none')
    set(H2,'LineStyle','-','LineWidth',2,'Color','k')
    plot(X1,zeros(1,length(X1)),'w','LineWidth',2,'LineStyle','--');
    %axis limits
    %     axis(AX(1),[0 length(X1) minY maxY])
    axis(AX(1),[0 length(X1) -500 500]);
    set(AX(1),'FontSize',16)
    set(AX(2),'FontSize',16,'YColor','k')
    axis(AX(2),[0 length(X1) 0 1])
    title({'Qualitatively forward, but quantitatively reversible. The group contribution ';
           '\Delta_{r}G^{\primem} estimates span the zero line.'},'FontSize',16)
    set(get(AX(1),'Ylabel'),'String','\Delta_{r}G^{\primem}  (blue)    or     \Delta_{r}G^{\prime} (red) (kJ/mol)')
    set(get(AX(2),'Ylabel'),'String','P(\Delta_{r}G^{\primem}<0)')
    set(get(AX(1),'Ylabel'),'FontSize',16)
    set(get(AX(2),'Ylabel'),'FontSize',16)
    xlabel('Reactions, sorted by P(\Delta_{r}G^{\primem}<0)');
end

changedDirections = directions;
% old code
%         %reactions that are qualitatively forward but quantitatively reverse
%         X=1:nnz(changedDirections.forwardReversible);
%         Y=(model.dGt0Min(changedDirections.forwardReversible)+model.dGt0Max(changedDirections.forwardReversible))/2;
%         L=Y-model.dGt0Min(changedDirections.forwardReversible);
%         U=model.dGt0Max(changedDirections.forwardReversible)-Y;
%         [tmp,xi]=sort(Y);
%         [tmp,xi]=sort(model.dGt0Max(changedDirections.forwardReversible)-model.dGt0Min(changedDirections.forwardReversible));
%         errorbar(X,Y(xi),L(xi),U(xi),'b.')
%         title('All reactions that are qualitatively forward, yet quantitatively reversible');
%         ylabel('dGt0 (kJ/mol)')
%         xlabel('Reactions, sorted by uncertainty in dGt0')
%
%         X=1:nnz(forwardReversibleKeq);
%         Y=(model.dGt0Min(forwardReversibleKeq)+model.dGt0Max(forwardReversibleKeq))/2;
%         L=Y-model.dGt0Min(forwardReversibleKeq);
%         U=model.dGt0Max(forwardReversibleKeq)-Y;
%         [tmp,xi]=sort(Y);
%         figure;
%         errorbar(X,Y(xi),L(xi),U(xi),'b.')
%         title('Reactions qualitiatively forward, but quantitatively reversible by Alberty');
%         ylabel('dGt0 (kJ/mol)')
%         xlabel('Reactions, sorted by mean dGt0')
% %     %GC forward prediction of reactions qualitiatively forward
% %     X1=1:nnz(dGfGCforwardForwardBool);
% %     Y=(model.dGt0Min(dGfGCforwardForwardBool)+model.dGt0Max(dGfGCforwardForwardBool))/2;
% %     L=Y-model.dGt0Min(dGfGCforwardForwardBool);
% %     U=model.dGt0Max(dGfGCforwardForwardBool)-Y;
% %     [tmp,xi]=sort(Y);
% %     figure;
% %     hold on
% %     errorbar(X1,Y(xi),L(xi),U(xi),'b.')
% %     legend('forwardForward')
% %     title('GC forward prediction of reactions qualitiatively forward')
% %     ylabel('dGt0 (kJ/mol)')
% %     xlabel('Reactions, sorted by mean dGt0')
% %
% %     [tmp,xi]=sort(abs(model.dGt0Max(dGfGCforwardForwardBool)-model.dGt0Min(dGfGCforwardForwardBool)));
% %     figure;
% %     hold on
% %     errorbar(X1,Y(xi),L(xi),U(xi),'b.')
% %     legend('forwardForward')
% %     title('GC forward prediction of reactions qualitiatively forward')
% %     ylabel('dGt0 (kJ/mol)')
% %     xlabel('Reactions, sorted by uncertainty in dGt0')
% %
% %     figure
% %     plot(Y(xi),L(xi),'.g')
% %     title('GC forward prediction of reactions qualitiatively forward')
% %     ylabel('dGt0 Uncertainty (kJ/mol)')
% %     xlabel('Reaction dGt0 (kJ/mol), sorted by uncertainty in dGt0')
%
%     %GC reversible prediction of reactions qualitiatively forward by dGrt
%     %dGr0t
%     X1=1:nnz(directions.forwardReversible);
%     Y=(model.dGt0Min(directions.forwardReversible)+model.dGt0Max(directions.forwardReversible))/2;
%     L=Y-model.dGt0Min(directions.forwardReversible);
%     U=model.dGt0Max(directions.forwardReversible)-Y;
% %     [tmp,xi]=sort(Y);
% %     figure
% %     errorbar(X1,Y(xi),L(xi),U(xi),'r.')
% %     title('All GC reversible prediction of reactions qualitiatively forward')
% %     legend('forwardReversible')
% %     ylabel('dGt0 (kJ/mol)')
% %     xlabel('Reactions, sorted by mean dGt0')
%     %dGrt
%     Y2=(model.dGtMin(directions.forwardReversible)+model.dGtMax(directions.forwardReversible))/2;
%     L2=Y2-model.dGtMin(directions.forwardReversible);
%     U2=model.dGtMax(directions.forwardReversible)-Y2;
%     [tmp,xi]=sort(abs(model.dGt0Max(directions.forwardReversible)-model.dGt0Min(directions.forwardReversible)));
%     figure
%     hold on;
%     %dGrt errorbar
%     hE=errorbar(X1,Y2(xi),L2(xi),U2(xi),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','r');
%     % adjust error bar width
% %     hE_c=get(hE, 'Children');
% %     errorbarXData= get(hE_c(2), 'XData');
%     errorbarXData= get(hE, 'XData');
%     errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0.2;
%     errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0.2;
%     %set(hE_c(2), 'XData', errorbarXData);
%     set(hE, 'XData', errorbarXData);
%     %dGrt0 errorbar on top and inside dGrt
%     hE2=errorbar(X1,Y(xi),L(xi),U(xi),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','b');
%     % adjust error bar width
%     %hE_c=get(hE2, 'Children');
%     %errorbarXData= get(hE_c(2), 'XData');
%     errorbarXData= get(hE2, 'XData');
%
%     errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0.2;
%     errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0.2;
%     %set(hE_c(2), 'XData', errorbarXData);
%     set(hE2, 'XData', errorbarXData);
%     %mean dGrt0
%     plot(X1,Y(xi),'.g','LineStyle','none');
%     %zero line
%     plot(X1,zeros(1,length(X1)),'k');
%     title('GC reversible prediction of reactions qualitiatively forward');
%     ylabel('dGt0 (kJ/mol)');
%     xlabel('Reactions, sorted by uncertainty in dGt0');
%
%     %%%%%%%%
%     %qualitatively forward reactions that are quantitatively reversible by concentration alone
%     %dGrt0
%     Y=(model.dGt0Min(dGfGCforwardReversibleBool_byConc)+model.dGt0Max(dGfGCforwardReversibleBool_byConc))/2;
%     L=Y-model.dGt0Min(dGfGCforwardReversibleBool_byConc);
%     U=model.dGt0Max(dGfGCforwardReversibleBool_byConc)-Y;
%     %dGrt
%     Y2=(model.dGtMin(dGfGCforwardReversibleBool_byConc)+model.dGtMax(dGfGCforwardReversibleBool_byConc))/2;
%     L2=Y2-model.dGtMin(dGfGCforwardReversibleBool_byConc);
%     U2=model.dGtMax(dGfGCforwardReversibleBool_byConc)-Y2;
%     [tmp,xi]=sort(abs(model.dGt0Max(dGfGCforwardReversibleBool_byConc)-model.dGt0Min(dGfGCforwardReversibleBool_byConc)));
%
%     figure
%     hold on;
%     %dGrt errorbar underneath
%     X1=1:nnz(dGfGCforwardReversibleBool_byConc);
%     hE=errorbar(X1,Y2(xi),L2(xi),U2(xi),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','r');
%     % adjust error bar width
%     %hE_c=get(hE, 'Children');
%     %errorbarXData= get(hE_c(2), 'XData');
%     errorbarXData= get(hE, 'XData');
%
%     errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0.2;
%     errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0.2;
%     %set(hE_c(2), 'XData', errorbarXData);
%     set(hE, 'XData', errorbarXData);
%
%     %dGrt0 errorbar on top and inside dGrt
%     hE2=errorbar(X1,Y(xi),L(xi),U(xi),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','b');
%     % adjust error bar width
%     %hE_c=get(hE2, 'Children');
%     %errorbarXData= get(hE_c(2), 'XData');
%     errorbarXData= get(hE2, 'XData');
%
%     errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0.2;
%     errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0.2;
%     %set(hE_c(2), 'XData', errorbarXData);
%     set(hE2, 'XData', errorbarXData);
%
%
%     %mean dGrt0
% %     plot(X1,Y(xi),'Marker','none','LineStyle','none');
%     plot(X1,Y(xi),'.g','LineStyle','none');
%     %zero line
%     plot(X1,zeros(1,length(X1)),'k');
%     title('Qualitatively forward reactions that are GC quantitatively reversible by concentration');
%     ylabel('dGt0 or dGt(kJ/mol)');
%     xlabel('Reactions, sorted by uncertainty in dGt0');
%
%
%     %qualitatively forward reactions that are quantitatively reversible by
%     %the range of dGt0
%     fprintf('%i%s\n',nnz(dGfGCforwardReversibleBool_bydGt0),' qualitatively forward reactions that are gc quantitatively reversible by range of dGt0');
%     X1=1:nnz(dGfGCforwardReversibleBool_bydGt0);
%     %dGrt0
%     Y=(model.dGt0Min(dGfGCforwardReversibleBool_bydGt0)+model.dGt0Max(dGfGCforwardReversibleBool_bydGt0))/2;
%     L=Y-model.dGt0Min(dGfGCforwardReversibleBool_bydGt0);
%     U=model.dGt0Max(dGfGCforwardReversibleBool_bydGt0)-Y;
%     %dGrt
%     Y2=(model.dGtMin(dGfGCforwardReversibleBool_bydGt0)+model.dGtMax(dGfGCforwardReversibleBool_bydGt0))/2;
%     L2=Y2-model.dGtMin(dGfGCforwardReversibleBool_bydGt0);
%     U2=model.dGtMax(dGfGCforwardReversibleBool_bydGt0)-Y2;
%     [tmp,xi]=sort(abs(model.dGt0Max(dGfGCforwardReversibleBool_bydGt0)-model.dGt0Min(dGfGCforwardReversibleBool_bydGt0)));
%
%     figure;
%     hold on;
%     %dGrt errorbar
%     hE=errorbar(X1,Y2(xi),L2(xi),U2(xi),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','r');
%     % adjust error bar width
%     %hE_c=get(hE, 'Children');
%     %errorbarXData= get(hE_c(2), 'XData');
%     errorbarXData= get(hE, 'XData');
%
%     errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0.2;
%     errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0.2;
%     %set(hE_c(2), 'XData', errorbarXData);
%     set(hE, 'XData', errorbarXData);
%
%     %dGrt0 errorbar on top and inside dGrt
%     hE2=errorbar(X1,Y(xi),L(xi),U(xi),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','b');
%     % adjust error bar width
%     %hE_c=get(hE2, 'Children');
%     %errorbarXData= get(hE_c(2), 'XData');
%     errorbarXData= get(hE2, 'XData');
%
%     errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0.2;
%     errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0.2;
%     %set(hE_c(2), 'XData', errorbarXData);
%     set(hE2, 'XData', errorbarXData);
%
%     %mean dGrt0
% %     plot(X1,Y(xi),'Marker','none','LineStyle','none');
%     plot(X1,Y(xi),'.g','LineStyle','none');
%     %zero line
%     plot(X1,zeros(1,length(X1)),'k');
%     title('Qualitatively forward reactions that are GC quantitatively reversible by range of dGt0');
%     ylabel('dGt0 or dGt(kJ/mol)');
%     xlabel('Reactions, sorted by uncertainty in dGt0');
%
%      %qualitatively forward reactions that are quantitatively reverseible by
%     %the range of dGt0
%     fprintf('%i%s\n',nnz(dGfGCforwardReversibleBool_bydGt0),' qualitatively forward reactions that are gc quantitatively reversible by range of dGt0');
%     X1=1:nnz(dGfGCforwardReversibleBool_bydGt0);
%     %dGrt0
%     Y=(model.dGt0Min(dGfGCforwardReversibleBool_bydGt0)+model.dGt0Max(dGfGCforwardReversibleBool_bydGt0))/2;
%     L=Y-model.dGt0Min(dGfGCforwardReversibleBool_bydGt0);
%     U=model.dGt0Max(dGfGCforwardReversibleBool_bydGt0)-Y;
%     %dGrt
%     Y2=(model.dGtMin(dGfGCforwardReversibleBool_bydGt0)+model.dGtMax(dGfGCforwardReversibleBool_bydGt0))/2;
%     L2=Y2-model.dGtMin(dGfGCforwardReversibleBool_bydGt0);
%     U2=model.dGtMax(dGfGCforwardReversibleBool_bydGt0)-Y2;
%     %find the amount of reactions with normal cumulative distribution
%     P = normcdf(0,Y,L);
%
%     %sort by fraction of uncertainty that is negative
% %     [tmp,xi]=sort(abs(model.dGt0Min(dGfGCforwardReversibleBool_bydGt0))./abs(model.dGt0Max(dGfGCforwardReversibleBool_bydGt0)-model.dGt0Min(dGfGCforwardReversibleBool_bydGt0)),'descend');
%     %     [tmp,xi]=sort(Y);
%     [tmp,xi]=sort(P,'descend');
%     %cuttoff for probabilities: must be reflective about 0.5;
%     PSortedeversibleBool=tmp<0.6 & tmp>0.4;
%     %reversible section bar in the background
%     figure1=figure;
%     hold on;
%     xbar=1:length(PSortedeversibleBool);
%     PSortedeversibleBar=ones(1,length(PSortedeversibleBool))*min(model.dGtMin(dGfGCforwardReversibleBool_bydGt0));
%     PSortedeversibleBar(PSortedeversibleBool)=max(model.dGtMax(dGfGCforwardReversibleBool_bydGt0));
%     bar_handle=bar(xbar,PSortedeversibleBar,1,'BaseValue',min(model.dGtMin(dGfGCforwardReversibleBool_bydGt0)),'EdgeColor',[0.86 0.86 0.86]);
%     %dGrt errorbar
%     hE=errorbar(X1,Y2(xi),L2(xi),U2(xi),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','r');
%     % adjust error bar width
%     %hE_c=get(hE, 'Children');
%     %errorbarXData= get(hE_c(2), 'XData');
%     errorbarXData= get(hE, 'XData');
%
%     errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0.2;
%     errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0.2;
%     %set(hE_c(2), 'XData', errorbarXData);
%     set(hE, 'XData', errorbarXData);
%
%     %dGrt0 errorbar on top and inside dGrt
%     hE2=errorbar(X1,Y(xi),L(xi),U(xi),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','b');
%     % adjust error bar width
%     %hE_c=get(hE2, 'Children');
%     %errorbarXData= get(hE_c(2), 'XData');
%     errorbarXData= get(hE2, 'XData');
%
%     errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0.2;
%     errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0.2;
%     %set(hE_c(2), 'XData', errorbarXData);
%     set(hE2, 'XData', errorbarXData);
%
%
%     %mean dGrt0
% %     plot(X1,Y(xi),'Marker','none','LineStyle','none');
%     plot(X1,Y(xi),'.g','LineStyle','none');
%     %zero line
%     %cumulative probability that reaction is really forward, assuming a
%     %normal distribution about the mean dGt0
%     [AX,H1,H2]=plotyy(X1,zeros(1,length(X1)),X1,P(xi));
%     set(AX(2),'YTick',[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1]);
%     set(H1,'LineStyle','none')
%     set(H2,'LineStyle','--','Color','k')
%     plot(X1,zeros(1,length(X1)),'k');
%     axis(AX(1),[0 length(X1) min(model.dGtMin(dGfGCforwardReversibleBool_bydGt0)) max(model.dGtMax(dGfGCforwardReversibleBool_bydGt0))])
%     axis(AX(2),[0 length(X1) 0 1])
%     title('Qualitatively forward reactions that are GC quantitatively reversible by range of dGt0');
%     ylabel('dGt0 or dGt(kJ/mol)');
%     xlabel('Reactions, sorted by cumulative probability of being forward');
%
%     %qualitatively forward reactions that are quantitatively reverse by
%     %dGt0
%     fprintf('%i%s\n',nnz(dGfGCforwardReversibleBool_byConcRHS),' qualitatively forward reactions that are GC quantitatively forward by dGf0');
%     X1=1:nnz(dGfGCforwardReversibleBool_byConcRHS);
%     %dGrt0
%     Y=(model.dGt0Min(dGfGCforwardReversibleBool_byConcRHS)+model.dGt0Max(dGfGCforwardReversibleBool_byConcRHS))/2;
%     L=Y-model.dGt0Min(dGfGCforwardReversibleBool_byConcRHS);
%     U=model.dGt0Max(dGfGCforwardReversibleBool_byConcRHS)-Y;
%     %dGrt
%     Y2=(model.dGtMin(dGfGCforwardReversibleBool_byConcRHS)+model.dGtMax(dGfGCforwardReversibleBool_byConcRHS))/2;
%     L2=Y2-model.dGtMin(dGfGCforwardReversibleBool_byConcRHS);
%     U2=model.dGtMax(dGfGCforwardReversibleBool_byConcRHS)-Y2;
%     [tmp,xi]=sort(abs(model.dGt0Max(dGfGCforwardReversibleBool_byConcRHS)-model.dGt0Min(dGfGCforwardReversibleBool_byConcRHS)));
%
%     figure
%     hold on;
%     %dGrt errorbar
%     hE=errorbar(X1,Y2(xi),L2(xi),U2(xi),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','r');
%     % adjust error bar width
%     %hE_c=get(hE, 'Children');
%     %errorbarXData= get(hE_c(2), 'XData');
%     errorbarXData= get(hE, 'XData');
%
%     errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0.2;
%     errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0.2;
%     %set(hE_c(2), 'XData', errorbarXData);
%     set(hE, 'XData', errorbarXData);
%
%     %dGrt0 errorbar on top and inside dGrt
%     hE2=errorbar(X1,Y(xi),L(xi),U(xi),'LineStyle','none','LineWidth',2,'DisplayName','forwardReversible','Color','b');
%     % adjust error bar width
%     %hE_c=get(hE2, 'Children');
%     %errorbarXData= get(hE_c(2), 'XData');
%     errorbarXData= get(hE2, 'XData');
%
%     errorbarXData(4:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(7:9:end) = errorbarXData(1:9:end) - 0.2;
%     errorbarXData(5:9:end) = errorbarXData(1:9:end) + 0.2;
%     errorbarXData(8:9:end) = errorbarXData(1:9:end) + 0.2;
%     %set(hE_c(2), 'XData', errorbarXData);
%     set(hE2, 'XData', errorbarXData);
%     %mean dGrt0
% %     plot(X1,Y(xi),'Marker','none','LineStyle','none');
%     plot(X1,Y(xi),'.g','LineStyle','none');
%     %zero line
%     plot(X1,zeros(1,length(X1)),'k');
%     title('Qualitatively forward reactions that are GC quantitatively forward by dGf0');
%     ylabel('dGt0 or dGt(kJ/mol)');
%     xlabel('Reactions, sorted by uncertainty in dGt0');

author: opencobra

repo: cobratoolbox


text:
% Q = vbfa(D, Y, W_module, X_module, noise_module, ...)
%
% Variational Bayesian (VB) factor analysis (FA) learning algorithm with
% changeable modules for the latent variables.

function Q = vbfa(D, Y, W_module, X_module, noise_module, varargin)

[M,N] = size(Y);

options = struct('maxiter', 100, ...
                 'update_x', 1, ...
                 'update_w', 1, ...
                 'update_noise', 1, ...
                 'rotate', 1, ...
                 'rotation_checkgrad', false, ...
                 'rotation_show', false, ...
                 'rotation_maxiter', 10, ...
                 'debug', false, ...
                 'autosave', false,...
                 'autosavefile', 'vbfa_autosave');

[options, errmsg] = argparse(options, varargin{:});
error(errmsg);

% Initialize
disp('Initializing variables..')
[Q.v_W,Q.W,Q.CovW] = W_module.initialize();
[Q.v_X,Q.X,Q.CovX] = X_module.initialize();
[Q.Tau] = noise_module.initialize();
Q.W_module = W_module;
Q.X_module = X_module;
Q.noise_module = noise_module;

% Observed values
Obs = ~isnan(Y);
MN_obs = sum(Obs(:));

% Replace missing values with zeros for computational reasons
Y(~Obs) = 0;

%%
%% VB learning

loglikelihood_old = -inf;
Q.loglikelihood = nan(options.maxiter,1);

Q.Y = Y;

KL_W = -inf;
KL_X = -inf;
KL_Tau = -inf;

%f = zeros(4,options.maxiter);

disp('Starting the VB iteration..')
for ind=1:options.maxiter

  startitercpu = cputime;

  %
  % Update variables
  %

  if index_selected(ind, options.update_x)
    if index_selected(ind, options.debug)
      disp('Update X')
    end
    [Q.v_X,Q.X,Q.CovX,KL_X] = X_module.update(ind, Y, Obs, Q.v_W, Q.W, Q.CovW, Q.Tau);
  end

  if index_selected(ind, options.update_w)
    if options.debug, disp('Update W'), end
    if index_selected(ind, options.debug)
      disp('Update W')
    end
    [Q.v_W,Q.W,Q.CovW,KL_W] = W_module.update(ind, Y', Obs', Q.v_X, Q.X, Q.CovX, Q.Tau');
  end

  % Compute squared errors <(y_mn-w_m*x_n)^2>
  % (used by noise update and loglikelihood)
  E2 = zeros(size(Y));
  Yh = Q.W'*Q.X;
  E2 = Y.*Y;
  E2 = E2 - 2*Y.*Yh;
  if ndims(Q.CovW) == 2 && ndims(Q.CovX) == 2
    E2 = E2 + Yh.^2 + Q.W'.^2*Q.CovX + Q.CovW'*Q.X.^2 + Q.CovW'*Q.CovX;
  elseif ndims(Q.CovW) == 3 && ndims(Q.CovX) == 3
    xx = bsxfun(@times, reshape(Q.X,[D,1,N]), reshape(Q.X,[1,D,N])) + Q.CovX;
    xx = reshape(xx, [D*D,N]);
    ww = bsxfun(@times, reshape(Q.W,[D,1,M]), reshape(Q.W,[1,D,M])) + Q.CovW;
    ww = reshape(ww, [D*D,M]);
    E2 = E2 + ww'*xx;
  else
    % TODO: Optimize this..
    warning('Optimize this..');
    for m=1:M
      for n=1:N
        if Obs(m,n)
          if ndims(Q.CovW) == 2
            ww = Q.W(:,m)*Q.W(:,m)' + diag(Q.CovW(:,m));
          else
            ww = Q.W(:,m)*Q.W(:,m)' + Q.CovW(:,:,m);
          end
          if ndims(Q.CovX) == 2
            xx = Q.X(:,n)*Q.X(:,n)' + diag(Q.CovX(:,n));
          else
            xx = Q.X(:,n)*Q.X(:,n)' + Q.CovX(:,:,n);
          end
          %WX_WX = WX_WX + traceprod(ww, xx);
          E2(m,n) = E2(m,n) + traceprod(ww, xx);
        end
      end
    end
  end
  E2(~Obs) = 0;

  if index_selected(ind, options.update_noise)
    if index_selected(ind, options.debug)
      disp('Update Tau')
    end

    [Q.Tau,LogTau,KL_Tau] = noise_module.update(ind, E2, Obs);
%    [Tau,lowerbound] = noise_module.update(Y, Obs, v_W, W, CovW, v_X, X, CovX);
  end

  %
  % Rotate
  %

  if index_selected(ind, options.rotate)

    % TODO: You could optimize the hyperparameters at the same time?
    disp('Rotating..')
    A = eye(D);

    if index_selected(ind, options.rotation_checkgrad)
      mycheckgrad(@rotation_cost, A(:) + 0.5*randn(D^2,1), 1e-3, W_module, ...
                  X_module);
    end
    A = minimize(A(:), @rotation_cost, options.rotation_maxiter, W_module, ...
                 X_module);
    A = reshape(A, [D D]);
    if index_selected(ind, options.rotation_show)
      A
    end
    [Q.W, Q.CovW] = W_module.rotate(A);
    [Q.X, Q.CovX] = X_module.rotate(inv(A)');

  end

  %
  % Evaluate VB lower bound
  %

  % Likelihood part: <log p(Y|...)>
  logpdf_y = gaussian_logpdf(Q.Tau(Obs)'*E2(Obs), ...
                             0, ...
                             0, ...
                             -sum(LogTau(Obs)), ...
                             MN_obs);

  % Lower bound
  Q.loglikelihood(ind) = logpdf_y - KL_W - KL_X - KL_Tau;

% $$$   if ~isreal(Q.loglikelihood(ind))
% $$$     KL_W
% $$$     KL_X
% $$$     KL_Tau
% $$$     error('Loglikelihood imaginary')
% $$$   end

  if Q.loglikelihood(ind) < loglikelihood_old
    warning(sprintf('Loglikelihood lower bound decreased relatively %e!', ...
                    (loglikelihood_old - Q.loglikelihood(ind)) / loglikelihood_old));
  end

  fprintf('Iteration step %d: loglikelihood=%e (%.2f seconds)\n', ind, ...
          Q.loglikelihood(ind), cputime()-startitercpu);

% $$$   f(1,ind) = KL_W;
% $$$   f(2,ind) = KL_X;
% $$$   f(3,ind) = KL_Tau;
% $$$   f(4,ind) = loglikelihood(ind);

  loglikelihood_old = Q.loglikelihood(ind);

  if index_selected(ind, options.autosave)
    fprintf('Saving results to %s...', options.autosavefile);
    save(options.autosavefile, '-struct', 'Q');
    fprintf(' done.\n');
  end

end

% $$$ figure
% $$$ plot(f');

% $$$ if nargout >= 1
% $$$
% $$$   Q.W_module = W_module;
% $$$   Q.W = W;
% $$$   Q.CovW = CovW;
% $$$
% $$$   Q.X_module = X_module;
% $$$   Q.X = X;
% $$$   Q.CovX = CovX;
% $$$
% $$$   Q.noise_module = noise_module;
% $$$   Q.Tau = Tau;
% $$$
% $$$   Q.loglikelihood = loglikelihood;
% $$$
% $$$ end

%function rotate(W_module, X_module)
%end

function [c, dc] = rotation_cost(a, W_module, X_module)
N = sqrt(length(a));
A = reshape(a, N, N);
[U,S,V] = svd(A);
invS = diag(1./diag(S));
invAt = U*invS*V';
[c_w, dc_w] = W_module.rotation_cost(A, U, S, V);
[c_x, dc_x] = X_module.rotation_cost(invAt, U, invS, V);
dc_x = -invAt*dc_x'*invAt;
c = c_w + c_x;
dc = dc_w(:) + dc_x(:);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [ y2, m2, d2, f2 ] = ymdf_inc_islamic ( y1, m1, d1, f1, days )

%*****************************************************************************80
%
%% YMDF_INC_ISLAMIC increments an Islamic YMDF date by DAYS days.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 April 2103
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1, M1, D1, real F1,
%    the YMDF date.
%
%    Input, real DAYS, the number of days to advance the date.
%
%    Output, integer Y2, M2, D2, real F2,
%    the incremented YMDF date.
%

%
%  Copy the parameters.
%
  y2 = y1;
  m2 = m1;
  d2 = d1;
  f2 = f1 + days;
%
%  Check the parameters.
%
  [ y2, m2, d2, f2, ierror ] = ymdf_check_islamic ( y2, m2, d2, f2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Demo script for comparing Matlab SVDS with LMSVD

clc
% problem size
%Asize = input('[m n] = ');
Asize = [1 1]*2000;
m = Asize(1);
n = Asize(2);

nR = 5; pcent = 5;
Ranks = (1:nR)*ceil(pcent/100*min(m,n)/nR);
ferr = @(u,v) norm(u(:)-v(:))/norm(v(:));

% generate random A
tau = rand;
base = 1 + 0.1*tau;
d = max(base.^(0:-1:-n+1),eps);
A = randn(m,n)*sparse(1:n,1:n,d);
A = A(randperm(m),randperm(n));

% set options
opts.tol = 1e-8;
opts.maxit = 150;

T1 = zeros(1,nR); E1 = T1;
T2 = T1; E2 = E1;

for j = 1:nR

    r = Ranks(j); k = r + 10;
    fprintf('\n size: [m n r] = [%i %i %i]\n',m,n,r)

    % Matlab svds
    tic; [U1,S1,V1] = svds(A,r,'L',opts);
    t1 = toc; e1 = ferr(U1*S1*V1',A); sv1 = diag(S1);
    fprintf(' svds: res = %16.12e, t = %5.3e\n',e1,t1);
    T1(j) = t1; E1(j) = e1;

    % lmsvd
    tic; [U2,S2,V2,Out] = lmsvd(A,r,opts);
    t2 = toc; e2 = ferr(U2*S2*V2',A); sv2 = diag(S2);
    fprintf('lmsvd: res = %16.12e, t = %5.3e\n',e2,t2);
    T2(j) = t2; E2(j) = e2;

end
fprintf('\nSummary: avg-err = %.2e. time ratio: %.2f\n\n',...
    mean(abs(E1-E2)),sum(T1)/sum(T2));

%% plot
close all; scrsz = get(0,'ScreenSize');
figure('Position',[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2])
subplot(121)
svk = Out.svk; bar(svk/svk(1),.01); hold on
plot(1:r,svk(1:r)/svk(1),'bo','MarkerSize',3,'linewidth',2);
%plot(r+1:k,svk(r+1:k)/svk(1),'bd','MarkerSize',3,'linewidth',2);
hold off
t_str = sprintf('Singular Values');
title(t_str,'fontsize',14);
xlabel('Number of sv''s','fontsize',14);
set(gca,'fontsize',12);
axis square

subplot(122)
semilogy(Ranks,T1,'rs:',Ranks,T2,'bd:',...
    'MarkerSize',6,'linewidth',2)
t_str = sprintf('Matrix: %i x %i',m,n);
title(t_str,'fontsize',14);
xlabel('Number of sv''s','fontsize',14);
ylabel('Time','fontsize',14);
legend('svds','lmsvd','location','Best');
set(gca,'fontsize',12); grid on
axis square
shg

author: indigits

repo: sparse-plex


text:
function obj = loadobj(s)
% LOADOBJ  Initialize a DagNN object from a structure.
%   OBJ = LOADOBJ(S) initializes a DagNN objet from the structure
%   S. It is the opposite of S = OBJ.SAVEOBJ().
%   If S is a string, initializes the DagNN object with data
%   from a mat-file S. Otherwise, if S is an instance of `dagnn.DagNN`,
%   returns S.

% Copyright (C) 2015 Karel Lenc and Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

if ischar(s) s = load(s); end
if isstruct(s)
  assert(isfield(s, 'layers'), 'Invalid model.');
  if ~isstruct(s.layers)
    warning('The model appears to be `simplenn` model. Using `fromSimpleNN` instead.');
    obj = dagnn.DagNN.fromSimpleNN(s);
    return;
  end
  obj = dagnn.DagNN() ;
  for l = 1:numel(s.layers)
    constr = str2func(s.layers(l).type) ;
    block = constr() ;
    block.load(struct(s.layers(l).block)) ;
    obj.addLayer(...
      s.layers(l).name, ...
      block, ...
      s.layers(l).inputs, ...
      s.layers(l).outputs, ...
      s.layers(l).params) ;
  end
  if isfield(s, 'params')
    for f = setdiff(fieldnames(s.params)','name')
      f = char(f) ;
      for i = 1:numel(s.params)
        p = obj.getParamIndex(s.params(i).name) ;
        obj.params(p).(f) = s.params(i).(f) ;
      end
    end
  end
  if isfield(s, 'vars')
    for f = setdiff(fieldnames(s.vars)','name')
      f = char(f) ;
      for i = 1:numel(s.vars)
        p = obj.getVarIndex(s.vars(i).name) ;
        obj.vars(p).(f) = s.vars(i).(f) ;
      end
    end
  end
  for f = setdiff(fieldnames(s)', {'vars','params','layers'})
    f = char(f) ;
    obj.(f) = s.(f) ;
  end
elseif isa(s, 'dagnn.DagNN')
  obj = s ;
else
  error('Unknown data type %s for `loadobj`.', class(s));
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function [ know, x ] = p01_sol ( m, know )

%*****************************************************************************80
%
%% P01_SOL returns the solution for problem 01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input/output, integer KNOW.
%    On input, KNOW is 0, or the index of the previously returned solution.
%    On output, KNOW is 0 if there are no more solutions, or it is the
%    index of the next solution.
%
%    Output, real X(M,1), the solution, if known.
%
  if ( know == 0 )
    know = 1;
    x(1:m,1) = 1.0;
  else
    know = 0;
    x(1:m,1) = 0.0;
  end

  return
end


text:
function varargout = Imagesegmentation(varargin)
%% Written by ABIOYE SAMSON D.
% Email Address: samsung4christ@yahoo.com
% Phone No: +2348068724645
%%
% IMAGESEGMENTATION M-file for Imagesegmentation.fig
%      IMAGESEGMENTATION, by itself, creates a new IMAGESEGMENTATION or raises the existing
%      singleton*.
%
%      H = IMAGESEGMENTATION returns the handle to a new IMAGESEGMENTATION or the handle to
%      the existing singleton*.
%
%      IMAGESEGMENTATION('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in IMAGESEGMENTATION.M with the given input arguments.
%
%      IMAGESEGMENTATION('Property','Value',...) creates a new IMAGESEGMENTATION or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before Imagesegmentation_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to Imagesegmentation_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

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

% Last Modified by GUIDE v2.5 27-Sep-2010 21:29:28

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

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

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

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

% Update handles structure
guidata(hObject, handles);

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

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

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
[filename path] = uigetfile('*.bmp;*.jpg;*.tif;','Open input image');
if filename~=0
    file = [path filename];
    orignalimg = imread(file);
    axes(handles.axes1);
    imshow(real(orignalimg));
else
    warndlg('Input image must be selected.','Warning');
end

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
j=getImage(handles.axes1);
jr=imresize(j,[512 512]);
jg=rgb2gray(jr);
h=histeq(jg);
b=qtdecomp(h,.1);
s=full(b);
blocks=repmat(double(0),size(s));
for dim=[512 256 128 64 32 16 8 4 2 1];
numblocks=length(find(s==dim));
if (numblocks > 0)
values = repmat(double(1),[dim dim numblocks]);
values(2:dim,2:dim,:)=0;
blocks=qtsetblk(blocks,s,dim,values);
end
end
blocks(end,1:end)=1;
blocks(1:end,end)=1;
f=(blocks());
axes(handles.axes2);
imshow(f)

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

% --- Executes on button press in pushbutton4.
function pushbutton4_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
exit = questdlg('Are you sure you want to Exit now','Exit','Yes','No','');
switch exit
    case 'Yes'
        delete(gcf)
    case 'No'
        return
end

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

author: Sable

repo: mcbench-benchmarks


text:
function [ labels ] = kernelkmeans_classifier(XGrid,X,centroid,eigens,kernel,kpar)
%KMEANS_CLASSIFIER K-means classifier
%
%  input ------------------------------------------------------------------
%
%       o X          : (N x D), data to classifiy
%
%       o centroids  : (L x D), L centroids
%
%       o dist       : string, distance metric
%
%   output-----------------------------------------------------------------
%
%       o labels     : (N x 1), class labels
%
%

KnGrid         = gram(XGrid,X,kernel,kpar(1),kpar(2));
KnGrid         = KnGrid * eigens;
HGrid          = KnGrid ./ repmat(sqrt(sum(KnGrid.^2, 2)), 1, size(eigens,2));

D           = ml_distfunc(HGrid, centroid, 'sqeuclidean');
[~, labels] = min(D, [], 2);

end

author: epfl-lasa


text:
function [h,g,a,info] = wfilt_oddevenb(N)
%WFILT_ODDEVENB  Kingsbury's symmetric odd filters
%
%   Usage: [h,g,a] = wfilt_oddevenb(N);
%
%   `[h,g,a]=wfilt_oddevenb(N)` with $N \in {1}$ returns Kingsbury's
%   odd filters.
%
%   Examples:
%   ---------
%   :::
%     figure(1);
%     wfiltinfo('ana:oddevenb1');
%
%     figure(2);
%     wfiltinfo('syn:oddevenb1');
%
%   References: king02

% AUTHOR: Zdenek Prusa

info.istight = 0;

a = [2;2];

switch(N)
 case 1
    % Example 1. from the reference. Symmetric near-orthogonal
    garr = [
             0              0
             0             -7.062639508928571e-05
             0              0
            -0.0017578125   1.341901506696429e-03
             0             -1.883370535714286e-03
             0.022265625   -7.156808035714285e-03
            -0.046875       2.385602678571428e-02
            -0.0482421875   5.564313616071428e-02
             0.2968750     -5.168805803571428e-02
             0.55546875    -2.997576032366072e-01
             0.2968750      5.594308035714286e-01
            -0.0482421875  -2.997576032366072e-01
            -0.046875      -5.168805803571428e-02
             0.022265625    5.564313616071428e-02
             0              2.385602678571428e-02
            -0.0017578125  -7.156808035714285e-03
             0             -1.883370535714286e-03
             0              1.341901506696429e-03
             0              0
             0             -7.062639508928571e-05
    ];

    % This scaling is not in the reference paper, but it is here to be
    % consistent
    garr = garr*sqrt(2);
    %garr = setnorm(garr,'energy');

    offset = -10;

  otherwise
        error('%s: No such filters.',upper(mfilename));
end

    %garr = [garr(:,3:4),garr(:,1:2)];
    modrange = (-1).^((0:size(garr,1)-1) + offset+1).';
    modrange2 = (-1).^((0:size(garr,1)-1) + offset).';

    harr =       [garr(:,2).*modrange2,...
                  garr(:,1).*modrange,...
                  ];

% In the biorthogonal case, the filters do not get time reversed
garr = flipud(garr);

htmp=mat2cell(harr,size(harr,1),ones(1,size(harr,2)));
h = cellfun(@(hEl)struct('h',hEl,'offset',offset),htmp(1:2),...
                   'UniformOutput',0);

gtmp=mat2cell(garr,size(garr,1),ones(1,size(garr,2)));

g = cellfun(@(gEl)struct('h',gEl,'offset',offset),gtmp(1:2),...
                   'UniformOutput',0);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat