text:
function [prob,opts] = buildOpti(varargin)
%Build an OPTI object
%
%   Called By OPTI Constructor

%   Copyright (C) 2011-2012 Jonathan Currie (IPL)

% Build Settings
settings.QUAD_EIG_MAX_SIZE = 1e8;    %Maximum size (rows*cols) for a quadratic matrix before skipping eig() for non-convex check

%API Change to make optiprob redundant (1/5/12)
[prob,opts] = exPrbOpts(varargin{:});

%Check and correct problem size and setup errors
[prob,opts] = checkOpti(prob,opts,settings);

function [prob,opts] = checkOpti(prob,opts,settings)
%Check the objective, constraints and options when building an OPTI object

%Get Warning Level
if(strcmpi(opts.warnings,'all'))
    warn = 2;
elseif(strcmpi(opts.warnings,'critical'))
    warn = 1;
else
    warn = 0;
end

%Numerical Differences Structure
numdif = struct('grad',0,'hess',0,'jac',0);
%AMPL Interface Structure
ampl = struct('path',prob.path,'useASL',0,'writesol',0);
%Misc Structure
misc = struct('Horig',[],'funorig',[],'ydata_orig',[],'xdata_orig',[],'weighting_orig',[],'forceSNLE',false,'fitFun',[],'fitGrad',[]);

%Sizes structure
siz = struct('ndec',[],'ncon',[],'nrow',[],'nineq',[],'neq',[],'nbnds',[],'nbin',[],'nint',[],'nqc',[],'nsdcone',[],'nnlrow',[],'nnlineq',[],'nnleq',[],'nsos',[],'nsos1',[],'nsos2',[]);

%--Check Objective Function--%
%Check we at least have an opti problem
if(isempty(prob.f) && isempty(prob.fun) && isempty(prob.ode) && (isempty(prob.sdcone) || ~isstruct(prob.sdcone) || ~isfield(prob.sdcone,'b')))
    %see if we have system of linear equations
    if(~isempty(prob.A) && ~isempty(prob.b) && isempty(prob.nlcon))
        [prob,opts] = buildSLE(prob,opts);
        prob.sizes = siz;
        prob.sizes.ndec = size(prob.A,1);
        prob.ampl = ampl;
        return;
    %see if we have system of nonlinear equations (to be solved as an NLP)
    elseif(~isempty(prob.nlcon) && isempty(prob.fun))
        misc.forceSNLE = true;
    else
    	error('You have not supplied the objective function in a readable form - or it is empty!');
    end
end

%Save original nonlinear objective
if(~isempty(prob.fun))
    misc.funorig = prob.fun;
end

%Transpose f if required
if(size(prob.f,2) > 1)
    prob.f = prob.f';
end
if(~isempty(prob.x0))
    [r,c] = size(prob.x0);
    if(r > 1 && c > 1)
        error('OPTI does not currently solve matrix problems. Please use reshape() within your objective / constraints, and pass x0 as a vector to OPTI to solve these problems');
    end
    if(c > 1)
        prob.x0 = prob.x0';
    end
    if(issparse(prob.x0))
        if(warn > 1), optiwarn('opti:sparsex0','The initial guess (x0) must be a dense vector'); end
        prob.x0 = full(prob.x0);
    end
end

%Check correct sizes for QP
if(~isempty(prob.H) && ~isa(prob.H,'function_handle'))
    if(isempty(prob.f))
        prob.f = zeros(size(prob.H,1),1);
    else
        [n,m] = size(prob.H);
        p = length(prob.f);
        if(n ~= p || m ~= p)
            error('The sizes of the QP H and f matrices don''t correspond!');
        end
    end
    misc.Horig = prob.H; %save for nnz later in display
end

%Get ndec if not NL
if(~isempty(prob.f) && ~isa(prob.f,'function_handle'))
    siz.ndec = length(prob.f);
elseif(~isempty(prob.H) && ~isa(prob.H,'function_handle'))
    siz.ndec = size(prob.H,1);
elseif(~isempty(prob.sdcone) && isstruct(prob.sdcone) && isfield(prob.sdcone,'b'))
    siz.ndec = length(prob.sdcone.b); %sedumi b
end

%Check x0 if we have ndec
if(~isempty(siz.ndec) && ~isempty(prob.x0))
    if(siz.ndec  ~= length(prob.x0))
        error('The supplied x0 is not the correct length, expected %d x 1',siz.ndec);
    end
end

%Check for NL grad
if(~isempty(prob.fun) && isempty(prob.f))
    %Get Solver Info
    info = optiSolverInfo(opts.solver,[],[],opts);
    %Use finite difference if requires derivative
    if(info.der1)
        if(warn > 1), optiwarn('opti:mkljac','OPTI will use MKLJAC (Numerical Difference Algorithm) for the Objective Gradient Function'); end
        %Check for Data fitting problem
        if(~isempty(prob.ydata))
            if (nargin(prob.fun) == 2)
                prob.f = @(x) mklJac(@(x) prob.fun(x,prob.xdata),x);
            else
                prob.f = @(x) mklJac(prob.fun,x,length(prob.ydata));
            end
        %Normal objective, single row
        else
            prob.f = @(x) mklJac(prob.fun,x,1);
        end
        numdif.grad = 1;
    end
end

%Check for NL hess (never really approximate using numdiff)
if(~isempty(prob.fun) && isempty(prob.f))
    numdif.hess = 1;
end

%--Check Constraints--%
%Get ndec
if(isempty(siz.ndec))
    if(~isempty(prob.x0))
        siz.ndec = length(prob.x0);
    elseif(~isempty(prob.lb))
        siz.ndec = length(prob.lb);
    elseif(~isempty(prob.ub))
        siz.ndec = length(prob.ub);
    elseif(~isempty(prob.f) && isnumeric(prob.f))
        siz.ndec = length(prob.f);
    elseif(~isempty(prob.H) && isnumeric(prob.H))
        siz.ndec = size(prob.H,1);
    elseif(~isempty(prob.A))
        siz.ndec = size(prob.A,2);
    elseif(~isempty(prob.Aeq))
        siz.ndec = size(prob.Aeq,2);
    elseif(~isempty(prob.int))
        siz.ndec = length(prob.int);
    elseif(~isempty(prob.nljacstr))
        siz.ndec = size(prob.nljacstr(),2);
    elseif(~isempty(prob.Hstr))
        siz.ndec = size(prob.Hstr(),2);
    else
        siz.ndec = []; %Can't determine sizes! All NL
    end
end

%Build in equalities if missing
if(misc.forceSNLE && isempty(prob.cl) && isempty(prob.nlrhs))
    if(isempty(siz.ndec))
        error('Please supply the number of decision variables via x0, or via cl,cu = zeros() if not square, to solve a system of nonlinear equations');
    end
    if(~isempty(prob.x0) && length(prob.x0) == siz.ndec && all(~isnan(prob.x0)))
        neqs = length(prob.nlcon(prob.x0));
    else
        neqs = length(prob.nlcon(zeros(siz.ndec,1)));
    end
    prob.cl = zeros(neqs,1);
    prob.cu = zeros(neqs,1);
end

%If sdcone supplied as a structure, assume SeDuMi format and convert to OPTI form (if not using sedumi as the solver)
if(~isempty(prob.sdcone) && isstruct(prob.sdcone) && isfield(prob.sdcone,'b'))
    if(strcmpi(opts.solver,'auto'))
        opts.solver = 'sedumi'; %use sedumi to solve sedumi problems
    end
    if(~strcmpi(opts.solver,'sedumi'))
        prob = sedumi2opti(prob); %otherwise convert to opti format
    end
end

%Process Integer Constraints
prob = checkInt(prob,siz.ndec);

%Determine what constraints are used
[cA,cb,cAeq,cbeq,cru,crl,clb,cub,cqc,cqrl,cqru,csdp,cint,cnl,ccl,ccu,csos] = conUsed(prob);

%Auto-fill constraint pairs if we know enough information
if(cA && crl && ~cru), prob.ru = Inf(size(prob.rl)); cru = 1; end
if(cA && cru && ~crl), prob.rl = -Inf(size(prob.ru)); crl = 1; end
if(cnl && ccl && ~ccu), prob.cu = Inf(size(prob.cl)); ccu = 1; end
if(cnl && ccu && ~ccl), prob.cl = -Inf(size(prob.cu)); ccl = 1; end

%Check correct constraint pairs used
if(crl && (cb || cbeq))
    error('Currently you cannot supply both linear row based constraints (rl <= Ax <= ru) and linear inequality constraints (Ax <= b, Aeqx = beq)');
end
if(cb && ~cA)
    error('You must supply A and b for Inequality Constraints!')
end
if(crl && ~cA)
    error('You must supply A, rl and ru for Linear Constraints!')
end
if(crl && length(prob.rl) ~= length(prob.ru))
    error('Constraint vectors rl and ru are not the same size');
end
if(xor(cAeq,cbeq))
    error('You must supply Aeq and beq for Equality Constraints!');
end
if(cnl)
    if(xor(isempty(prob.cl),isempty(prob.cu)))
        error('You must supply both bounds (cl, cu) for Nonlinear Constraints!');
    end
    if(isempty(prob.cl) && isempty(prob.nlrhs))
        %If we know ndec, we can assume all constraints are <= 0, otherwise error
        if(~isempty(siz.ndec))
            %See if we have nle
            if(~isempty(prob.nle))
                if(warn)
                   optiwarn('opti:nlnorhs','You have not supplied the right hand side (nlrhs) or bounds (cl,cu) of the Nonlinear Constraints. OPTI will assume nlrhs = zeros(nnlcon,1).');
                end
                nnl = length(prob.nle);
                prob.nlrhs = zeros(nnl,1);
            else
                if(warn)
                   optiwarn('opti:nlempty','You have not supplied the right hand side (nlrhs), nonlinear equation type (nle) or bounds (cl,cu) of the Nonlinear Constraints. OPTI will assume all nonlinear constraints are <= 0.');
                end
                testC = prob.nlcon(zeros(siz.ndec,1)); nnl = length(testC);
                prob.nlrhs = zeros(nnl,1);
                prob.nle = -1*ones(nnl,1);
            end
        else
            error('You must supply the Right Hand Side (nlrhs) OR constraint bounds (cl, cu) for Nonlinear Constraints! If you wish to default to all <= 0, please supply ndec.');
        end
    end
    if(~isempty(prob.cl) && ~isempty(prob.nlrhs))
        error('Currently you cannot supply both nonlinear row based constraints (cl <= nlcon(x) <= cu) and nonlinear mixed constraints (nlcon, nlrhs, nle)');
    end
end
if(~isempty(prob.sos) && xor(isempty(prob.sos.index),isempty(prob.sos.weight)))
    error('You must supply the SOS type, indices and weights!');
end

%See if we are constrained or not
if(~cA && ~cAeq && ~clb && ~cub && ~cqc && ~csdp && ~cint && ~cnl && ~csos)
    prob.iscon = 0;
    checkLP(prob);
else
    prob.iscon = 1;
end

%Get Constraint Sizes
if(crl)
    eq = prob.rl == prob.ru;
    neq = sum(eq);
    nineq = sum(~isinf(prob.rl))+sum(~isinf(prob.ru))-2*neq;
    nrow = length(prob.rl);
else
    if(cb)
        nineq = length(prob.b);
    elseif(isstruct(prob.sdcone) && isfield(prob.sdcone,'K') && isfield(prob.sdcone.K,'l'))
        nineq = prob.sdcone.K.l;
    else
        nineq = 0;
    end
    if(cbeq); neq = length(prob.beq); else neq = 0; end
    nrow = 0;
end
if(clb); nbnds = length(find(~isinf(prob.lb) > 0)); else nbnds = 0; end
if(cub); nbnds = nbnds + length(find(~isinf(prob.ub) > 0)); end
if(cint)
    nbin = length(find(prob.int.str == 'B'));
    nint = length(find(prob.int.str == 'I'));
else
    nbin = 0; nint = 0;
end
if(cnl)
    if(~isempty(prob.nlrhs))
        nnl = length(prob.nlrhs);
    else
        nnl = length(prob.cl);
    end
    nnlrow = length(prob.cl);
else
    nnl = 0;
    nnlrow = 0;
end
if(cqc)
    if(iscell(prob.Q)), nqc = length(prob.Q); else nqc = 1; end
else
    nqc = 0;
end
if(csdp)
    if(iscell(prob.sdcone))
        nsdcone = length(prob.sdcone);
    elseif(isstruct(prob.sdcone) && isfield(prob.sdcone,'K') && isfield(prob.sdcone.K,'s'))
        nsdcone = length(prob.sdcone.K.s);
    else
        nsdcone = 1;
    end
else
    nsdcone = 0;
end
if(csos)
    nsos = length(prob.sos.type);
    [r,c] = size(prob.sos.type);
    if(r > c)
        t = str2num(prob.sos.type); %#ok<*ST2NM>
    else
        t = str2num(prob.sos.type');
    end
    nsos1 = sum(t==1);
    nsos2 = sum(t==2);
else
    nsos = 0;
    nsos1 = 0;
    nsos2 = 0;
end

%Collect number of constraints
ncon = nineq + neq + nbnds + nqc + nsdcone + nbin + nint + nnl + nsos;

%Transpose as neccesary
if(cb && size(prob.b,2) > 1)
    prob.b = prob.b';
end
if(cbeq && size(prob.beq,2) > 1)
    prob.beq = prob.beq';
end
if(crl && size(prob.rl,2) > 1)
    prob.rl = prob.rl';
end
if(cru && size(prob.ru,2) > 1)
    prob.ru = prob.ru';
end
if(clb && size(prob.lb,2) > 1)
    prob.lb = prob.lb';
end
if(cub && size(prob.ub,2) > 1)
    prob.ub = prob.ub';
end
if(cnl && size(prob.nlrhs,2) > 1)
    prob.nlrhs = prob.nlrhs';
end
if(cnl && size(prob.cl,2) > 1)
    prob.cl = prob.cl';
end
if(cnl && size(prob.cu,2) > 1)
    prob.cu = prob.cu';
end
%Check and transpose quadratic constraints
if(cqc)
    %Convert l, r to numeric if cells
    if(iscell(prob.l))
        if(size(prob.l,1) > 1 && size(prob.l,2) > 1)
            error('Quadratic constraint l must be a column cell array');
        end
        n = length(prob.l);
        l = zeros(siz.ndec,n);
        for i = 1:n
            li = prob.l{i};
            if(size(li,1) > 1 && size(li,2) > 1)
                error('Quadratic constraint l in cell %d must be a column vector!',i);
            end
            if(size(li,2) > 1)
                li = li';
            end
            l(:,i) = li;
        end
        prob.l = l;
    end
    prob.l = chkdouble(prob.l,'prob.l',warn > 1);
    if(iscell(prob.qrl))
        if(size(prob.qrl,1) > 1 && size(prob.qrl,2) > 1)
            error('Quadratic constraint qrl must be a column cell array');
        end
        n = length(prob.qrl);
        qrl = zeros(n,1);
        for i = 1:n
            qrli = prob.qrl{i};
            if(size(qrli,1) > 1 || size(qrli,2) > 1)
                error('Quadratic constraint qrl in cell %d must be a column scalar!',i);
            end
            qrl(i) = qrli;
        end
        prob.qrl = qrl;
    end
    if(iscell(prob.qru))
        if(size(prob.qru,1) > 1 && size(prob.qru,2) > 1)
            error('Quadratic constraint qru must be a column cell array');
        end
        n = length(prob.qru);
        qru = zeros(n,1);
        for i = 1:n
            qrui = prob.qru{i};
            if(size(qrui,1) > 1 || size(qrui,2) > 1)
                error('Quadratic constraint qru in cell %d must be a column scalar!',i);
            end
            qru(i) = qrui;
        end
        prob.qru = qru;
    end
    prob.qrl = chkdouble(prob.qrl,'prob.qrl',warn > 1);
    prob.qru = chkdouble(prob.qru,'prob.qru',warn > 1);
    if(iscell(prob.Q))
        if(size(prob.Q,2) > 1)
            prob.Q = prob.Q';
        end
        if(numel(prob.Q)==1)
            prob.Q = prob.Q{1};
        end
    else
        if(size(prob.l,2) > 1)
            prob.l = prob.l';
        end
    end
    if(size(prob.qrl,2) > 1)
        prob.qrl = prob.qrl';
    end
    if(size(prob.qru,2) > 1)
        prob.qru = prob.qru';
    end
end
%Auto fill missing quadratic bounds
if(cqc && cqrl && ~cqru), prob.qru = Inf(size(prob.qrl)); cqru = 1; end
if(cqc && cqru && ~cqrl), prob.qrl = -Inf(size(prob.qru)); cqrl = 1; end
if(cqrl && ~cqc)
    error('You must supply Q, l, qrl and qru for Quadratic Constraints!')
end
if(cqrl && length(prob.qrl) ~= length(prob.qru))
    error('Constraint vectors qrl and qru are not the same size');
end

%Check and transpose SOS constraints
if(csos)
    if(size(prob.sos.type,2) > 1)
        prob.sos.type = prob.sos.type';
    end
    if(iscell(prob.sos.index))
        if(size(prob.sos.index,2) > 1)
            prob.sos.index = prob.sos.index';
        end
        for i = 1:length(prob.sos.index)
            if(size(prob.sos.index{i},2) > 1)
                prob.sos.index{i} = prob.sos.index{i}';
            end
        end
    else
        if(size(prob.sos.index,2) > 1)
            prob.sos.index = prob.sos.index';
        end
    end
    if(iscell(prob.sos.weight))
        if(size(prob.sos.weight,2) > 1)
            prob.sos.weight = prob.sos.weight';
        end
        for i = 1:length(prob.sos.weight)
            if(size(prob.sos.weight{i},2) > 1)
                prob.sos.weight{i} = prob.sos.weight{i}';
            end
        end
    else
        if(size(prob.sos.weight,2) > 1)
            prob.sos.weight = prob.sos.weight';
        end
    end
end

%Fill in nle if empty (assumes all <=)
if(cnl && isempty(prob.cl) && isempty(prob.nle))
    if(warn)
        optiwarn('opti:build','Nonlinear Constraint Bounds (cl, cu) and Type (nle) are Empty - Assuming all Nonlinear Constraints are nlcon(x) <= nrhs');
    end
    prob.nle = -1*ones(nnl,1);
end

%Check for NL jac
if(cnl && isempty(prob.nljac))
    %Get Solver Info
    info = optiSolverInfo(opts.solver,[],[],opts);
    %Use finite difference if requires derivative
    if(info.der1)
        if(warn > 1),optiwarn('opti:mkljac','OPTI will use MKLJAC (Numerical Difference Algorithm) for the Constraint Jacobian Function'); end
        numdif.jac = 1;
        prob.nljac = @(x) mklJac(prob.nlcon,x,nnl);
    end
end

%Check Sizes if possible
if(cA && (siz.ndec ~= size(prob.A,2)))
    error('Constraint A matrix is the wrong size! Expected %d x %d (nineq x ndec)',nineq,siz.ndec);
end
if(cAeq && (siz.ndec ~= size(prob.Aeq,2)))
    error('Constraint Aeq matrix is the wrong size! Expected %d x %d (neq x ndec)',neq,siz.ndec);
end
if(clb && (siz.ndec ~= length(prob.lb)))
    error('Incorrect size of lb! Expected %d x 1 (ndec x 1)',siz.ndec);
end
if(cub && (siz.ndec ~= length(prob.ub)))
    error('Incorrect size of ub! Expected %d x 1 (ndec x 1)',siz.ndec);
end
if(crl)
    if(length(prob.rl) ~= size(prob.A,1))
        error('Constraint A matrix is the wrong size! Expected %d x %d (nlincon x ndec)',nineq+neq,siz.ndec);
    end
    if(length(prob.rl) ~= length(prob.ru))
        error('Linear constraints bounds rl and ru are not the same length!');
    end
else
    if(cA && isempty(prob.b))
        error('You must supply both A and b for inequality constraints');
    end
    if(cA && (length(prob.b) ~= size(prob.A,1)))
        error('Constraint A matrix is the wrong size! Expected %d x %d (nineq x ndec)',nineq,siz.ndec);
    end
    if(cAeq && isempty(prob.beq))
        error('You must supply both Aeq and beq for inequality constraints');
    end
    if(cAeq && (length(prob.beq) ~= size(prob.Aeq,1)))
        error('Constraint Aeq matrix is the wrong size! Expected %d x %d (neq x ndec)',neq,siz.ndec);
    end
end
if(cnl && (length(prob.nlrhs) ~= length(prob.nle)))
    error('Nonlinear RHS and e vectors are not the same length!');
end
if(ccl && ccu && (length(prob.cl) ~= length(prob.cu)))
    error('Nonlinear cl and cu vectors are not the same length!');
end
if(cqc)
    if(iscell(prob.Q))
        if(length(prob.Q) ~= size(prob.l,2) || siz.ndec ~= size(prob.l,1))
            error('Quadratic Constraint l matrix is the wrong size! Expected %d x %d (ndec x nqc)',siz.ndec,length(prob.Q));
        end
        if(length(prob.Q) ~= length(prob.qrl))
            error('Quadratic Constraint qrl vector is the wrong size! Expected %d x %d (nqc x 1)',length(prob.Q),1);
        end
    else
        if(1 ~= size(prob.l,2) || siz.ndec ~= size(prob.l,1))
            error('Quadratic Constraint l vector is the wrong size! Expected %d x %d (ndec x nqc)',siz.ndec,1);
        end
        if(1 ~= length(prob.qrl))
            error('Quadratic Constraint qrl vector is the wrong size! Expected %d x %d (nqc x 1)',1,1);
        end
    end
end
if(csos)
    ns = length(prob.sos.type);
    if(ns > 1)
        if(~iscell(prob.sos.index))
            error('SOS Index argument is not a cell! Expected a cell array of %d x 1 (nsos x 1)',ns);
        elseif(length(prob.sos.index) ~= ns)
            error('SOS Index cell array is the wrong size! Expected %d x 1 (nsos x 1)',ns);
        end
        if(~iscell(prob.sos.weight))
            error('SOS Weights argument is not a cell! Expected a cell array of %d x 1 (nsos x 1)',ns);
        elseif(length(prob.sos.weight) ~= ns)
            error('SOS Weights cell array is the wrong size! Expected %d x 1 (nsos x 1)',ns);
        end
        for i = 1:ns
            if(length(prob.sos.index{i}) ~= length(prob.sos.weight{i}))
                error('SOS Constraint %d Index and Weight vectors are not the same length!',i);
            end
            if(any(prob.sos.index{i} < 1))
                error('SOS Constraint %d Index vector must contain indices greater than 0 (MATLAB Indexing)!',i);
            end
        end
    else
        if(iscell(prob.sos.index) || iscell(prob.sos.weight))
            error('Only use cell arrays for multiple SOS constraints. Either sosind or soswt is a cell array - but there is only one SOS.');
        end
        if(length(prob.sos.index) ~= length(prob.sos.weight))
            error('SOS Index and Weight vectors are not the same length!');
        end
        if(any(prob.sos.index < 1))
            error('SOS Index vector must contain indices greater than 0 (MATLAB Indexing)!');
        end
    end
end

%-- Assign opti specified solver if present --%
if(isfield(prob,'solver') && ~isempty(prob.solver) && strcmpi(opts.solver,'auto'))
    opts.solver = prob.solver;
end
if(isfield(prob,'solver')), prob = rmfield(prob,'solver'); end

% Check Quadratic Constraints
QCisconvex = true;
if(cqc)
    %If we have any quadratic equalities OR double sided inequalities (assumed non-convex)
    %and the user has not specified a solver, default to SCIP (if available)
    if(any(~isinf(prob.qrl) & ~isinf(prob.qru)))
        if(strcmpi(opts.solver,'auto') && optiSolver('scip',0))
            opts.solver = 'scip';
        end
    end
    if(issparse(prob.l))
       prob.l = full(prob.l); %ensure dense
    end
    if(iscell(prob.Q))
       if(iscell(prob.l) || iscell(prob.qrl) || iscell(prob.qru))
           error('Only Quadratic Constraints Q may be a cell!');
       end
       r0 = size(prob.Q,1);
       c1 = size(prob.l,2);
       r2 = size(prob.qrl,1);
       r3 = size(prob.qru,1);
       if(r0 ~= c1 || r0 ~= r2 || r0 ~= r3)
           error('Quadratic Constraints Q + l + qrl + qru are not the same length!');
       end
       for i = 1:r0
           %Check Data Type
           prob.Q{i} = chkdouble(prob.Q{i},['prob.Q{',i,'}'],warn > 1);
           Q = prob.Q{i};
           l = prob.l(:,i);
           rl = prob.qrl(i);
           ru = prob.qru(i);
           %Check constraint
           QCisconvex = QCisconvex & chkQC(Q,l,rl,ru,i,siz.ndec,opts.solver,settings.QUAD_EIG_MAX_SIZE);
       end
    else
        %Check Data Type
        prob.Q = chkdouble(prob.Q,'prob.Q',warn > 1);
        %Check constraint
        QCisconvex = QCisconvex & chkQC(prob.Q,prob.l,prob.qrl,prob.qru,1,siz.ndec,opts.solver,settings.QUAD_EIG_MAX_SIZE);
    end
end

% Check Semidefinite Constraints
if(csdp)
    if(~isfield(siz,'ndec') || isempty(siz.ndec))
        error('In order to process semidefinite constraints you must supply ndec as well');
    end
    if(iscell(prob.sdcone))
        %Check Data type
        for i = 1:length(prob.sdcone)
            prob.sdcone{i} = chkdouble(prob.sdcone{i},sprintf('prob.sdcone{%d}',i),warn > 1);
            %Check and convert representation
            prob.sdcone{i} = chkSDCone(prob.sdcone{i},i,siz.ndec,warn);
        end
    elseif(isstruct(prob.sdcone)) %sedumi format
        %Correct structure names
        if(isfield(prob.sdcone,'A'))
            prob.sdcone.At = prob.sdcone.A;
            prob.sdcone = rmfield(prob.sdcone,'A');
        end
        if(isfield(prob.sdcone,'C'))
            prob.sdcone.c = prob.sdcone.C;
            prob.sdcone = rmfield(prob.sdcone,'C');
        end
        %Check data types
        prob.sdcone.At = chkdouble(prob.sdcone.At,'prob.sdcone.At',warn > 1);
        prob.sdcone.b = chkdouble(prob.sdcone.b,'prob.sdcone.b',warn > 1);
        prob.sdcone.c = chkdouble(prob.sdcone.c,'prob.sdcone.c',warn > 1);
        if(isfield(prob.sdcone,'K') && ~isstruct(prob.sdcone.K))
            error('SeDuMi prob.sdcone.K should be a structure!');
        end
    else
        %Check data type
        prob.sdcone = chkdouble(prob.sdcone,'prob.sdcone',warn > 1);
        %Check and convert representation
        prob.sdcone = chkSDCone(prob.sdcone,1,siz.ndec,warn);
    end
end

% Check Special Ordered Sets
if(csos)
    %Check types are correct
    types = str2num(prob.sos.type);
    if(any(types < 1 | types > 2))
        error('Only SOS types 1 and 2 are allowed!');
    end
    %Check Datatypes
    if(nsos > 1)
        %Check vectors
        for i = 1:nsos
            %Check Data Types
            prob.sos.index{i} = chkdouble(prob.sos.index{i},['prob.sos.index{',i,'}'],warn > 1);
            prob.sos.weight{i} = chkdouble(prob.sos.weight{i},['prob.sos.weight{',i,'}'],warn > 1);
        end
    else
        prob.sos.index = chkdouble(prob.sos.index,'prob.sos.index',warn > 1);
        prob.sos.weight = chkdouble(prob.sos.weight,'prob.sos.weight',warn > 1);
    end
end

%Check bounds direction
if(clb && cub)
    if(any(prob.ub < prob.lb))
        error('A lower bound is greater than its respective upper bound!');
    end
end
if(ccl && ccu)
    if(any(prob.cu < prob.cl))
        error('A nonlinear constraint lower bound is greater than its respective upper bound!');
    end
end
if(crl && cru)
    if(any(prob.ru < prob.rl))
        error('A linear constraint lower bound is greater than its respective upper bound!');
    end
end
if(cqrl && cqru)
    if(any(prob.qru < prob.qrl))
        error('A quadratic constraint lower bound is greater than its respective upper bound!');
    end
end

%Check for sparse bounds (causes problems when size doesn't match nz entries)
if(clb && issparse(prob.lb))
    error('The lower bound vector (lb) must be dense');
end
if(cub && issparse(prob.ub))
    error('The upper bound vector (ub) must be dense');
end
if(crl && issparse(prob.rl))
    error('The linear constraint lower bound vector (rl) must be dense');
end
if(cru && issparse(prob.ru))
    error('The linear constraint upper bound vector (ru) must be dense');
end
if(cb && issparse(prob.b))
    error('The linear inequality constraint RHS (b) must be dense');
end
if(cbeq && issparse(prob.beq))
    error('The linear equality constraint RHS (beq) must be dense');
end
if(ccl && issparse(prob.cl))
    error('The nonlinear constraint lower bound vector (cl) must be dense');
end
if(ccu && issparse(prob.cu))
    error('The nonlinear constraint upper bound vector (cu) must be dense');
end

%Check nonlinear constraints
if(cnl)
    if(~ccl && (any(prob.nle > 1) || any(prob.nle < -1)))
        error('Nonlinear constraint type (nle) must only contain -1, 0 or 1');
    end
end

%Fill in Sizes Structure
siz.ncon = ncon;
siz.nineq = nineq;
siz.neq = neq;
siz.nbnds = nbnds;
siz.nbin = nbin;
siz.nint = nint;
siz.nqc = nqc;
siz.nsdcone = nsdcone;
siz.nsos = nsos;
siz.nsos1 = nsos1;
siz.nsos2 = nsos2;
siz.nrow = nrow;
siz.nnlrow = nnlrow;
if(cnl)
    if(ccl)
        eq = prob.cl == prob.cu; nteq = ~eq;
        ile = isfinite(prob.cu) & nteq;
        ige = isfinite(prob.cl) & nteq;
        siz.nnlineq = sum(ile + ige);
        siz.nnleq = sum(eq);
        %If we have dual bounds, our number of con will be out
        if(any(ile & ige))
            siz.ncon = siz.ncon + sum(ile & ige);
        end
    else
        siz.nnlineq = sum(prob.nle ~= 0);
        siz.nnleq = sum(prob.nle == 0);
    end
else
    siz.nnlineq = 0;
    siz.nnleq = 0;
end
prob.sizes = siz;
%Fill in numerical differences structure
prob.numdif = numdif;

%-- Determine Problem Type --%
%Pre Check
if(isempty(prob.H) && isempty(prob.fun) && isempty(prob.ode) && ~misc.forceSNLE)
    prb = 'LP';
elseif(isempty(prob.ode) && isempty(prob.fun) && ~misc.forceSNLE)
    prb = 'QP';
elseif(isempty(prob.ode))
    prb = 'NLP';
else
    prb = 'DNLP';
end
%Check we don't have function handles for H, f
if(~strcmpi(prb,'NLP') && ~strcmp(prb,'DNLP'))
    if(isa(prob.H,'function_handle')) %QP
        error('When solving a Quadratic Problem the H matrix must be a double matrix - not a function handle');
    end
    if(isa(prob.f,'function_handle')) %LP or %QP
        error('When solving a Linear or Quadratic Problem the f vector must be a double vector - not a function handle');
    end
end
%Check we don't have a constant objective bias with nlp/nls
if(any(strcmpi(prb,{'DNLP','NLP'})) && isfield(prob,'objbias') && prob.objbias ~= 0)
    error('A constant objective bias is only supported with Linear and Quadratic Programs');
end
%Check for quadratic constraints
if(nqc)
    %Check if we have an LP
    if(strcmpi(prb,'LP'))
        prob.H = spalloc(length(prob.f),length(prob.f),0);
        prb = 'QP';
    %Otherwise if not a QP
    elseif(~strcmpi(prb,'QP'))
        error('Currently Quadratic Constraints are only supported with Quadratic Programs. Please reformulate your problem as a general nonlinear problem with nonlinear constraints');
    end
    prb = ['QC' prb];
end
if(nsdcone)
    if(~strcmpi(prb,'LP'))
        error('Currently only Semidefinite Constraints are supported with Linear Programs');
    end
    prb = 'SDP';
end
%Check for integer variables
if(any(prob.int.ind ~= 0) || nsos)
    if(all(prob.int.ind < 0) && strcmpi(prb,'LP')) %only allow BILPs
        prb = ['BI' prb];
    else
        prb = ['MI' prb];
    end
end
%Check if NLS
if((~isempty(prob.xdata) || ~isempty(prob.ydata)))
    if(strcmpi(prb,'NLP'))
        prb = 'NLS';
    elseif(strcmpi(prb,'DNLP'))
        prb = 'DNLS';
    else
        error('Currently only Nonlinear Least Squares Problems are Supported');
    end
end
%DNLPs not yet supported
if(strcmpi(prb,'DNLP'))
    error('Dynamic Nonlinear Programs are not yet supported');
end

%Check if UNO (scalar) or SNLE/SCNLE (vector)
if(strcmpi(prb,'NLP') && ~misc.forceSNLE)
    % Get x0 (or try and make one)
    if(isfield(prob,'x0') && ~isempty(prob.x0))
        x0 = prob.x0;
    elseif(prob.sizes.ndec)
        x0 = zeros(prob.sizes.ndec,1);
        if(warn > 1)
            optiwarn('OPTI:nox0','OPTI is testing your objective function but does not have an x0 to use - assuming zeros(n,1).\n Please supply x0 to opti/optiprob to avoid this warning.');
		end
    else
        error('OPTI cannot determine whether you are solving a UNO or SNLE. Please supply ndec or x0 to opti/optiprob to continue.');
    end
    try
        testObj = prob.fun(x0);
        l = length(testObj);
    catch ME
        fprintf(2,'OPTI detected an error running a test objective call. Please correct the error below and try again:\n\n');
        rethrow(ME); %so user can see it
    end
    % Ensure objective returns a double
    if (~isa(testObj, 'double'))
        error('OPTI:NotDouble', 'The objective function does not return a double precision value (%s instead)', class(testObj));
    end
    if(l > 1)
        %Check n = n
        if(l ~= length(x0))
            error('OPTI can only solve nonlinear equations which have n equations and n variables');
        end
		%If constrained, ensure only linearly or bound constrained
		if(prob.iscon)
			if(cqc || cnl || cint || csos)
				error('Only Linear Constraints and Bounds are supported with Nonlinear Equation Solving Problems');
			end
			prb = 'SCNLE';
		else
			prb = 'SNLE';
		end
        %Correct for numdiff size (defaults to scalar objective)
        if(numdif.grad)
            prob.f = @(x) mklJac(prob.fun,x,l);
        end
    elseif(~prob.iscon)
        prb = 'UNO';
    end
end

%Check we don't have nlcon with lp/qp/other
if(cnl && isempty(strfind(prb,'NLP')))
    error('Nonlinear Constraints are only supported with Nonlinear Programs. Please reformulate your problem as a general nonlinear problem with nonlinear constraints');
end

%Force SNLE/SCNLE type (based on missing fun to opti, but has nlcon)
if(misc.forceSNLE)
    if(cA || clb || cub || cint)
        prb = 'SCNLE';
    else
        prb = 'SNLE';
    end
    %Check neqs = nvars
%     if((prob.sizes.nnleq+prob.sizes.neq) > prob.sizes.ndec)
%         error('OPTI can only solve nonlinear equations which have n equations and n variables');
%     end
    %Check user hasn't supplied grad as well
    if(~isempty(prob.f))
        error(['OPTI is confused what problem you are trying to solve. When you don''t supply an objective and only nonlinear constraints, '...
               'OPTI assumes this is a SNLE. However you have also supplied an objective gradient..? Either supply the SNLE using ''fun'' and ''grad'''...
               ' or as a series of constraints using ''nlcon'' and ''nljac''']);
    end
    if(isempty(siz.ndec))
        error('OPTI requires the number of decisions variables to continue pre-processing. Please supply ndec or x0 to opti()');
    end
    %Check for sparse derivatives, if so, solve as NLP, otherwise, convert to normal OPTI SNLE problem
    isNormalSNLE = true; haveDeriv = false; isJacSparse = false;
    if(~isempty(prob.nljac))
        haveDeriv = true;
        %Assume if we have sparsity pattern then it is also sparse
        if(~isempty(prob.nljacstr))
            isNormalSNLE = false; isJacSparse = true;
        %Test if we have a sparse Jacobian
        else
            if(~isempty(prob.x0) && length(prob.x0) == siz.ndec && all(~isnan(prob.x0)))
                test_x0 = prob.x0;
            else
                test_x0 = zeros(siz.ndec,1);
            end
            if(issparse(prob.nljac(test_x0)))
                isNormalSNLE = false; isJacSparse = true;
            end
        end
    end
    %Check for NLP solver specified to solve problem
    info = optiSolverInfo(opts.solver);
    if(~strcmpi(opts.solver,'auto') && info.con.nlineq)
        isNormalSNLE = false; %assume nl constraints indicate NLP solver
    end
    %Normal problem
    if(isNormalSNLE)
        %Ensure only nonlinear equalities which = 0
        if(prob.sizes.nnlineq > 0)
            error('OPTI only solves SNLE problems with nonlinear EQUALITY constraints');
        end
        if((~isempty(prob.nlrhs) && any(prob.nlrhs ~= 0)) || (~isempty(prob.cu) && any(prob.cu ~= 0)))
            error('OPTI only solves SNLE problems with a right hand side that = 0');
        end
        %Transfer problem over
        prob.fun = prob.nlcon;
        prob.nlcon = []; prob.nlrhs = []; prob.nle = []; prob.cl = []; prob.cu = [];
        prob.sizes.ncon = prob.sizes.ncon - prob.sizes.nnleq;
        prob.sizes.nnleq = 0; prob.sizes.nnlrow = 0;
        if(haveDeriv)
            prob.f = prob.nljac; prob.nljac = []; prob.nljacstr = [];
        end
    else %convert to dummy nlp
        prob.fun = @(x) -1;
        prob.f = @(x) zeros(1,siz.ndec);
        %Change solver to IPOPT/Bonmin if auto
        if(strcmpi(opts.solver,'auto'))
            if(any(prob.int.str ~= 'C')) %not really allowed?
                opts.solver = 'bonmin';
            else
                opts.solver = 'ipopt';
            end
        else
            %Check user has chosen a sparse NLP solver
            info = optiSolverInfo(opts.solver);
            if(isJacSparse && ~info.sparse)
                error(['As you have supplied sparse derivatives OPTI has converted your SNLE/SCNLE problem into an NLP. However your selected solver (%s) '...
                       'does not support sparse problems. Please choose an NLP solver that supports sparsity, such as IPOPT.'],opts.solver);
            end
        end
    end
end
%Save misc structure
prob.misc = misc;

%Check for qc with NLP solver
if(cnl && cqc)
	error('Specifying separate quadratic and nonlinear constraints is not currently supported. Please specify all constraints as nonlinear constraints');
end

%Matlab MI check
if(strcmpi(opts.solver,'matlab'))
    % If the user has specified MATLAB as the solver, ensure it is available
    % for the problem type (optiSolver can't check this before here)
    if (~any(strcmpi(optiSolver(prb),'matlab')))
        solver = optiSolver(['best_' prb]);
        if (warn)
            optiwarn('opti:nomatlab','MATLAB''s Optimization Toolbox is not available.\nUsing %s instead.',solver);
        end
        opts = optiset(opts,'solver',solver);
    elseif(any(strcmpi(prb,{'MILP','MIQP','MIQCQP','SDP','MINLP'})))
        if(~strcmpi(prb,'milp') && isempty(which('intlinprog.m'))) %check for 2014a intlinprog
            solver = optiSolver(['best_' prb]);
            if(warn)
                optiwarn('opti:mi','MATLAB does not currently supply a Mixed Integer or Semidefinite Solver\nUsing %s instead.',solver);
            end
            opts = optiset(opts,'solver',solver);
        end
    end
end

%If user has specified a problem type manually, compare to auto one
if(isfield(prob,'probtype') && ~isempty(prob.probtype))
    if(~strcmpi(prob.probtype,prb))
        %Check user field exists
        ptypes = optiSolver('ptypes');
        switch(lower(prob.probtype))
            case ptypes
                if(warn > 1)
                    optiwarn('opti:ptype','You have specified the problem type as a %s but OPTI identified the problem type as a %s. OPTI''s decision will be overridden.',upper(prob.probtype),upper(prb));
                end
            otherwise
                error('Unknown problem type specified in optiprob: %s',upper(prob.probtype));
        end
    end
    prb = upper(prob.probtype);
end

%Save problem type
prob.type = prb;
prob.solvetype = prb; %what the problem will be solved as
if(strcmp(prob.Name,'OPTI Problem'))
    prob.Name = ['OPTI ' upper(prb) ' Problem'];
end

%Check correct number of nonlinear constraints (requires a call to nlcon but error otherwise is pretty hard to follow)
if(cnl && ~isempty(prob.nlcon))
    % Get x0 (or try and make one)
    if(isfield(prob,'x0') && ~isempty(prob.x0))
        x0 = prob.x0;
    else
        x0 = zeros(prob.sizes.ndec,1);
        if(warn > 1)
            optiwarn('OPTI:nox0','OPTI is testing your nonlinear constraint function(s) but does not have an x0 to use - assuming zeros(n,1).\n Please supply x0 to opti/optiprob to avoid this warning.');
        end
    end
    testNLCon = prob.nlcon(x0);
    nnlcon = length(testNLCon);
    % Ensure constraints return a double
    if (~isa(testNLCon, 'double'))
        error('OPTI:NotDouble', 'The nonlinear constraint function does not return a double precision value (%s instead)', class(testNLCon));
    end

    %Determine number of constraint equations from bounds (remember double sided constraints will cause issues otherwise)
    if(~isempty(prob.cl))
        tnnlcon = length(prob.cl);
    else
        tnnlcon = length(prob.nle);
    end
    if(nnlcon ~= tnnlcon)
        error('The number of elements returned from the nonlinear constraints does not match the number of expected constraints, based on length(nlrhs/nle/cl/cu). Expected %d, got %d.',tnnlcon,nnlcon)
    end
end

%-- Check (and Correct) Options based on Problem Type & Settings --%

%BILP with bounds
switch(prb)
    case 'BILP'
        if((~isempty(prob.lb) || ~isempty(prob.ub)) && warn > 1)
            optiwarn('opti:bilp','Bounds are not currently used in BILP problems');
        end
end

%LPSOLVE instead of LP_SOLVE
if(strcmpi(opts.solver,'lpsolve'))
    opts.solver = 'lp_solve';
end

%Check for lmder variants
if(any(strcmpi(opts.solver,{'lm_der','lm_dif','lmdif'})))
    opts.solver = 'lmder';
end

%CBC for LPs
if(strcmpi(prb,'LP') && strcmpi(opts.solver,'cbc'))
    opts.solver = 'clp';
    if(warn > 1)
        optiwarn('OPTI:LPCBC','CBC is for solving MILPs, using CLP instead.');
    end
end

%Check constraints on bounded and linearly constrained solvers
if((siz.nineq + siz.neq + siz.nnleq + siz.nnlineq) > 0)
    if(strcmpi(opts.solver,'lbfgsb'))
        if(warn > 1)
            optiwarn('OPTI:BoundedNLP','L-BFGS-B can only solve bounded problems - using IPOPT instead.')
        end
        opts.solver = 'ipopt';
    end
end
if((siz.nnleq + siz.nnlineq) > 0)
    if(strcmpi(opts.solver,'pswarm'))
        if(warn > 1)
            optiwarn('OPTI:LinconNLP','PSwarm can only solve linearly constrained problems - using IPOPT instead.')
        end
        if(misc.forceSNLE)
            error('PSwarm cannot solve Systems of Nonlinear Equations, please choose another solver');
        else
            opts.solver = 'ipopt';
        end
    end
end
if(((siz.neq + siz.nnleq) > 0) && ((siz.nint + siz.nbin) > 0)) %MI with Equalities
    if(strcmpi(opts.solver,'gmatlab'))
        error('The MATLAB Global Optimization GA Algorithm with Integer Constraints cannot solve Equality Constrained Problems - Please use another solver.');
    end
end

% %Bonmin for NLPs
% if((strcmpi(prb,'UNO') || strcmpi(prb,'NLP')) && strcmpi(opts.solver,'bonmin'))
%     opts.solver = 'ipopt';
%     if(warn > 1)
%         optiwarn('OPTI:NLPBONMIN','BONMIN is for solving MINLPs, using IPOPT instead.');
%     end
% end

%Unconstrained QP (-H\f)
if(strcmpi(prb,'QP') && ~prob.iscon)
    opts.solver = 'matlab';
    if(warn > 1)
        optiwarn('OPTI:UnconQP','You have supplied an Unconstrained QP. This will be solved using -H\\f.');
    end
end

%Use MKLTRNLS for DNLS
if(strcmpi(opts.solver,'auto') && strcmpi(prb,'DNLS'))
    opts.solver = 'mkltrnls';
end

%AUTO solver - replace with best for problem type
if(strcmpi(opts.solver,'auto'))
    try
        opts.solver = lower(optiSolver(['best_' prb]));
    catch ME
        %Convert QCQP problems to general nonlinear (for now)
        if(any(strcmpi(prb,{'QCQP','MIQP','MIQCQP'})))
            if(warn)
                optiwarn('OPTI:NoQCQP','OPTI cannot find a solver to solve a %s explicitly. It will converted to a (MI)NLP to be solved.',upper(prb));
            end
            %Do conversion then return (internally calls this function)
            [prob,opts] = QCQP2NLP(prob,opts);
            prob.solvetype = 'NLP';
            return;
        else
            rethrow(ME);
        end
    end
end

%Only certain solvers solve SOS problems
if(csos && ~any(strcmpi(opts.solver,{'CPLEX','LP_SOLVE','CBC','SCIP','BONMIN'})))
    oS = opts.solver;
    if(optiSolver('cplex',0)) %prefer cplex, but may not be available
        opts.solver = 'cplex';
    elseif(optiSolver('scip',0))
        opts.solver = 'scip';
    elseif(optiSolver('cbc',0))
        opts.solver = 'cbc';
    elseif(optiSolver('lp_solve',0))
        opts.solver = 'lp_solve';
    else
        error('You do not have a solver which can solve problems specified with SOS');
    end
    if(warn > 1)
        optiwarn('OPTI:NoSOS','%s does not solve MILPs with Special Ordered Sets. Using %s instead.',upper(oS),upper(opts.solver));
    end
end

%Check weighting is for NLS or DNLS only
if(~isempty(prob.weighting) && ~any(strcmpi(prb,{'NLS','DNLS'})))
    error('Data fitting weights are only supported for nonlinear least squares (NLS) and parameter estimation problems (DNLS)');
end

%Check data fitting and constrained equation solving problems
if(any(strcmpi(prb,{'NLS','DNLS','SCNLE'})))
    %Check ydata & weighting is a vector
    if(~isempty(prob.ydata))
        if(size(prob.ydata,1) > 1 && size(prob.ydata,2) > 1 && ~iscell(prob.ydata))
            prob.ydata = prob.ydata(:);
        end
        if(~isempty(prob.weighting) && size(prob.weighting,1) > 1 && size(prob.weighting,2) > 1)
            prob.weighting = prob.weighting(:);
        end
        %Transpose as required
        if(size(prob.ydata,2) > size(prob.ydata,1))
            prob.ydata = prob.ydata';
        end
        if(~isempty(prob.xdata) && size(prob.xdata,2) > size(prob.xdata,1))
            prob.xdata = prob.xdata';
        end
        if(~isempty(prob.weighting) && size(prob.weighting,2) > size(prob.weighting,1))
            prob.weighting = prob.weighting';
        end
        %TAKE THE SQRT of WEIGHTS (matches Stats toolbox + Curve Fitting toolbox)
        if(iscell(prob.weighting))
            for i = 1:length(prob.weighting)
                prob.weighting{i} = sqrt(prob.weighting{i});
            end
        else
            prob.weighting = sqrt(prob.weighting);
        end
    end

    if(~isempty(prob.weighting))
        %Check ydata length = weighting length
        if(length(prob.ydata) ~= length(prob.weighting))
            if(iscell(prob.ydata))
                error('OPTI expects the same number of cells for both ydata and weighting when solving data fitting problems');
            else
                error('OPTI expects one weight per measurement (ydata) when solving data fitting problems');
            end
        end
        %Check weights aren't inf or nan
        if(iscell(prob.weighting))
            for i = 1:length(prob.weighting)
                if(any(isnan(prob.weighting{i})) || any(isinf(prob.weighting{i})))
                    error('Data fitting weights cannot be Inf or NaN!');
                end
            end
        else
            if(any(isnan(prob.weighting)) || any(isinf(prob.weighting)))
                error('Data fitting weights cannot be Inf or NaN!');
            end
        end
    end

	%Check correct solver for constrained problems
	if(prob.iscon)
		oS = opts.solver;
		sInfo = optiSolverInfo(opts.solver);
		%Check for linearly constrained
		if((neq+nineq) > 0 && (sInfo.con.lineq + sInfo.con.leq) == 0)
			rub = optiSolver('levmar'); %#ok<NASGU> %only levmar solves linearly constrained NLS
            opts.solver = 'levmar';
			if(warn > 1)
				if(any(strcmpi(prb,{'NLS','DNLS'})))
					optiwarn('OPTI:LinConNLS','%s does not solve Linearly Constrained Nonlinear Least Squares Problems. Using %s instead.',upper(oS),upper(opts.solver));
				else
					optiwarn('OPTI:LinConNLE','%s does not solve Linearly Constrained Nonlinear Equation Solving Problems. Using %s instead.',upper(oS),upper(opts.solver));
				end
			end
		%Check for bounded problems
		elseif(nbnds && ~sInfo.con.bnd)
			rub = optiSolver('mkltrnls'); %#ok<NASGU>
			opts.solver = 'mkltrnls';
			if(warn > 1)
				if(any(strcmpi(prb,{'NLS','DNLS'})))
					optiwarn('OPTI:BndNLS','%s does not solve Bounded Nonlinear Least Squares Problems. Using %s instead.',upper(oS),upper(opts.solver));
				else
					optiwarn('OPTI:BndNLE','%s does not solve Bounded Nonlinear Equation Solving Problems. Using %s instead.',upper(oS),upper(opts.solver));
				end
			end
		end
    end
    %Check for dynamic parameter fitting problems
    if(~isempty(prob.ode))
        prob = DNLS2NLS(prob,opts);
    elseif(iscell(prob.ydata) || iscell(prob.xdata))
        error('Cell arrays for xdata/ydata/weighting are only supported for dynamic parameter estimation problems');
    else
        %Save fitting functions for use in statistics routines
        prob.misc.fitFun = prob.fun;
        if(numdif.grad == 1 && nargin(prob.fun)==2)
            prob.misc.fitGrad = @(x,xdata) mklJac(@(x) prob.fun(x,xdata),x);
        elseif(isempty(prob.f))
            %Try build gradient as above
            if(nargin(prob.fun)==2)
                prob.misc.fitGrad = @(x,xdata) mklJac(@(x) prob.fun(x,xdata),x);
            else
                prob.misc.fitGrad = @(x) mklJac(prob.fun,x);
            end
        else
            prob.misc.fitGrad = prob.f;
        end
    end
end

%Check xdata supplied for fun(x,xdata) problems (NOT DNLS)
if(strcmpi(prb,'NLS') && nargin(prob.fun) > 1 && isempty(prob.xdata))
    error('The supplied objective function requires more than 1 input, but xdata is empty! OPTI assumes fun(x,xdata) based on this configuration.');
end

%Binary constraints for OPTI, BONMIN, GMATLAB - use bounds
if(any(strcmpi(opts.solver,{'opti','bonmin','gmatlab'})) && any(prob.int.ind < 0))
    %lower bounds
    if(isempty(prob.lb))
        prob.lb = -Inf(siz.ndec,1);
        prob.lb(prob.int.ind < 0) = 0;
    else
        lb = prob.lb(prob.int.ind < 0);
        lb(lb < 1) = 0; %only change out of binary bounds ones (could be 1 <= x <= 1)
        prob.lb(prob.int.ind < 0) = lb;
    end
    %upper bounds
    if(isempty(prob.ub))
        prob.ub = Inf(siz.ndec,1);
        prob.ub(prob.int.ind < 0) = 1;
    else
        ub = prob.ub(prob.int.ind < 0);
        ub(ub > 0) = 1;
        prob.ub(prob.int.ind < 0) = ub;
    end
end

%Check sense, add minus if maximization [I assume this is ok for non-convex problems too?]
if(prob.sense < 0)
    switch(prb)
        case {'LP','MILP','BILP','SDP','MISDP'}
            prob.f = -prob.f;
        case {'QP','QCQP','MIQP','MIQCQP'}
            prob.H = -prob.H;
            prob.f = -prob.f;
        case {'UNO','NLS','NLP','MINLP'}
            prob.fun = @(x) -prob.fun(x);
            if(~isempty(prob.f))
                prob.f = @(x) -prob.f(x);
            end
            if(~isempty(prob.H))
               switch(nargin(prob.H))
                   case 1
                       prob.H = @(x) -prob.H(x); %with no lambda assume no constraints, minus is ok
                   case 3
                       prob.H = @(x,sigma,lambda) prob.H(x,-sigma,lambda); %with sigma we can invert the sign
                   otherwise
                       error('OPTI cannot modify the Hessian Callback in order to maximise the problem. Please manually reformulate your problem as a maximization one.');
               end
            end
        case 'DNLS'
            error('OPTI cannot modify the ''sense'' in a dynamic parameter fitting problem');
    end
end

%Check if the problem was read via AMPL
if(~isempty(ampl.path))
    %If we are solving via SCIP, skip MATLAB callbacks (whitebox solver)
    if(any(strcmp(opts.solver,{'scip'})))
        ampl.useASL = true;
        ampl.writesol = false; %done internally within scip
        %Close asl interface if still open
        if(asl('isopen')), asl('close'); end
    else
        ampl.useASL = false;
        ampl.writesol = true;
    end
else
    ampl.useASL = false;
    ampl.writesol = false;
end
%Assign ampl structure
prob.ampl = ampl;
%Assign empty save structure
prob.save = [];

%Check for non-convex QP/QCQP
if(~isempty(strfind(prb,'QP')))
    if(prod(size(prob.H)) < settings.QUAD_EIG_MAX_SIZE) %#ok<PSIZE>
        try
           e = eig(prob.H);
        catch %#ok<CTCH>
           e = eigs(prob.H);
        end
        %Check for non-convex QP
        if(any(e < 0))
            if(strcmpi(opts.solver,'cplex') && ~isfield(opts.solverOpts,'solutiontarget'))
               opts.solverOpts.solutiontarget.Cur = 2; %allow indefinite
            elseif(warn)
                optiwarn('OPTI:NonconvexQP','QP objective appears non-convex, you may not obtain a globally optimal solution');
            end
        end
    end
    %Check for non-convex QCs
    if(~QCisconvex)
        if(strcmpi(opts.solver,'cplex'))
           opts.solverOpts.solutiontarget.Cur = 3; %allow indefinite
        elseif(warn)
            optiwarn('OPTI:NonconvexQC','QCQP constraint(s) appear(s) non-convex, you may not obtain a globally optimal solution');
        end
    end
end

%Remove unused fields
if(isfield(prob,'probtype')), prob = rmfield(prob,'probtype'); end

%Run Derivative Checker If Requested
if(strcmpi(opts.derivCheck,'on'))
    checkDerivs(prob,warn);
end

function c = checkInt(c,ndec)
%Check and assign the Integer Constraints

if(isempty(c.int))
    c.int = struct('str',char('C'*ones(1,ndec)),'ind',zeros(1,ndec),'idx',[]);
    return;
end

if(~ischar(c.int) && ~isnumeric(c.int))
    error('The integer variable string must be a char array or double vector!');
end

if(ischar(c.int) && (length(c.int) ~= ndec))
    error('Integer Variable String is not the correct length! Expected 1 x %d',ndec);
end
if(isnumeric(c.int) && (length(c.int) > ndec))
    error('Integer Variable Index Vector is not the correct length! Expected MAX 1 x %d',ndec);
end

str = ones(1,ndec);
ind = zeros(1,ndec);

if(ischar(c.int))
    lc = lower(c.int);
    cidx = lc=='c';
    bidx = lc=='b';
    iidx = lc=='i';
    str(cidx) = 'C'; ind(cidx) = 0;
    str(bidx) = 'B'; ind(bidx) = -1;
    str(iidx) = 'I'; ind(iidx) = 1;
    idx = find(iidx);
    if(sum(cidx) + sum(bidx) + sum(iidx) ~= length(c.int))
        error('Unknown character(s) in Integer Variable String! OPTI recognises ''C'' continuous, ''I'' integer and ''B'' binary.');
    end
else
    if(any(c.int == 0))
        error('Indicies of integer variables must be >= 1')
    end
    if(any(c.int > ndec))
        error('An integer index is > ndec!');
    end
    str = char('C' * str); str(c.int) = 'I';
    ind(c.int) = 1;
    idx = c.int;
end

c.int = struct('str',char(str),'ind',ind,'idx',idx);

function [cA,cb,cAeq,cbeq,cru,crl,clb,cub,cqc,cqrl,cqru,csdp,cint,cnl,ccl,ccu,csos] = conUsed(c)
%Determine what constraints are not empty

if(isempty(c.A)); cA = 0; else cA = 1; end
if(isempty(c.b)); cb = 0; else cb = 1; end
if(isempty(c.Aeq)); cAeq = 0; else cAeq = 1; end
if(isempty(c.beq)); cbeq = 0; else cbeq = 1; end
if(isempty(c.ru)); cru = 0; else cru = 1; end
if(isempty(c.rl)); crl = 0; else crl = 1; end
if(isempty(c.lb)); clb = 0; else clb = 1; end
if(isempty(c.ub)); cub = 0; else cub = 1; end
if(isempty(c.Q)); cqc = 0; else cqc = 1; end
if(isempty(c.qrl)); cqrl = 0; else cqrl = 1; end
if(isempty(c.qru)); cqru = 0; else cqru = 1; end
if(isempty(c.sdcone)); csdp = 0; else csdp = 1; end
if(any(c.int.str ~= 'C')); cint = 1; else cint = 0; end
if(isempty(c.nlcon)); cnl = 0; else cnl = 1; end
if(isempty(c.cl)); ccl = 0; else ccl = 1; end
if(isempty(c.cu)); ccu = 0; else ccu = 1; end
if(isempty(c.sos) || ~isfield(c.sos,'type') || isempty(c.sos.type)); csos = 0; else csos = 1; end

function checkLP(prob)
%Check if we are solving an unconstrained LP
if(isempty(prob.H) && isempty(prob.fun) && isempty(prob.ode) && ~prob.iscon)
    error('When solving an LP the problem must be constrained!');
end

function [prob,opts] = buildSLE(prob,opts)
%Build Problem for System of Linear Equations

prob.type = 'SLE';
prob.solvetype = 'SLE'; %what the problem will be solved as

%Check we have A + b + correct sizes
if(isempty(prob.A) || isempty(prob.b))
    error('You must supply both A + b to solve a SLE');
end
if(size(prob.A,1) ~= length(prob.b))
    error('Equation RHS vector b is the wrong size! Expected %d x 1',size(prob.A,1));
end

%AUTO solver - replace with best for problem type
if(strcmpi(opts.solver,'auto'))
    opts.solver = lower(optiSolver(['best_' prob.type]));
end

%MLDIVIDE solver
if(strcmpi(opts.solver,'mldivide'))
    opts.solver = 'matlab';
end

%no constraints
prob.iscon = 0;

function [prob,opts] = exPrbOpts(varargin)
%Extract the problem and options from the supplied arguments to opti

switch(nargin)
    case 0
        optiprob;
        error('You cannot create an empty OPTI object');
    case 1 %opti(prob) or opti(optiObj)
        if(isempty(varargin{1}))
            error('You cannot create an empty OPTI object');
        elseif(isstruct(varargin{1}))
            prob = optiprob(varargin{1});
            opts = optiset(); %default opts
        elseif (isa(varargin{1},'optim.problemdef.OptimizationProblem'))
            [prob,opts] = convertMatlabOptimProbToOptiProb(varargin{1});
        elseif(isa(varargin{1},'opti'))
            error('opti(optiObj) is not implemented');
        elseif(ischar(varargin{1}))
            if(~isempty(strfind(varargin{1},'.'))) %filename
                prob = optiRead(varargin{1});
                opts = optiset();
            else
                error('Unknown constructor calling form with argument ''%s''.\n- If you intended to read a mathematical file, you must supply the file extension.',varargin{1});
            end
        else
            error('Unknown calling form for OPTI!');
        end
    case 2 %opti(prob,opts) OR opti('field',value) OR opti(optiObj,opts) OR opti('file',opts)
        if(isstruct(varargin{1}) && isstruct(varargin{2}))
            prob = optiprob(varargin{1});
            opts = optiset(varargin{2});
        elseif (isa(varargin{1},'optim.problemdef.OptimizationProblem'))
            [prob,opts] = convertMatlabOptimProbToOptiProb(varargin{1},varargin{2});
        elseif(ischar(varargin{1}))
            if(~isempty(strfind(varargin{1},'.'))) %filename
                if(isstruct(varargin{2}))
                    prob = optiRead(varargin{1});
                    opts = optiset(varargin{2});
                else
                    error('Unknown calling form, expected to read a mathematical file with the second arg specified as an optiset structure');
                end
            else
                prob = optiprob(varargin{:});
                opts = optiset();
            end
        elseif(isa(varargin{1},'opti'))
            if(isstruct(varargin{2}) && isfield(varargin{2},'solver'))
                prob = getProb(varargin{1});
                opts = optiset(varargin{2});
            else
                error('opti(optiObj,arg) is not implemented unless arg is an options structure');
            end
        else
            error('Unknown calling form for OPTI!');
        end

    otherwise
        if(isa(varargin{1},'opti'))
            oprob = getProb(varargin{1});
        else
            oprob = [];
        end

%         try
            if(~isempty(oprob))
                prob = optiprob(oprob,varargin{2:end});
            else
                prob = optiprob(varargin{:});
            end
%         catch ME
%             throw(ME);
%         end
        if(~isempty(prob.opts))
            opts = prob.opts;
        else
            opts = optiset(); %default
        end
        %If solver is specified in prob and is auto in opts, set it
        if(strcmpi(opts.solver,'auto') && isfield(prob,'solver') && ~isempty(prob.solver))
            opts.solver = prob.solver;
            prob = rmfield(prob,'solver'); %don't keep a copy
        end
end
if(isfield(prob,'opts'))
    prob = rmfield(prob,'opts'); %don't keep a copy
end

function var = chkdouble(var,name,warn)
%Check if a variable is a double, if not, correct
if(~isa(var,'double'))
    if(warn)
        optiwarn('opti:notDouble','Variable %s was not a double, correcting: [double(%s)]',name,name);
    end
    var = double(var);
end

function isconvex = chkQC(Q,l,rl,ru,i,ndec,solver,QUAD_EIG_MAX_SIZE)
%Check a quadratic constraint for solver suitability
isconvex = true;
if((ndec ~= size(Q,1)) || (ndec ~= size(Q,2)))
    error('Constraint Q matrix in cell %d is the wrong size! Expected %d x %d',i,ndec,ndec);
end
if(ndec ~= size(l,1))
   error('Constraint l vector in column %d is the wrong size! Expected %d x 1',i,ndec);
end
if(~isscalar(rl))
   error('Constraint qrl in column %d is the wrong size! Expected 1 x 1 ',i);
end
if(~isscalar(ru))
   error('Constraint qru in column %d is the wrong size! Expected 1 x 1 ',i);
end
if(~any(strcmpi(solver,{'auto','scip','baron','ipopt','bonmin'}))) %allow scip or NLP solver to solve any QC constraint form
    %Check for equality
    if(rl == ru)
        error('Quadratic Constraint %d is an equality constraint. Please specify a global solver such as SCIP to solve this problem.\nAlternatively you may try IPOPT/BONMIN for a local solution.',i);
    end
    %If within 'sensible' size limit for checking convexity
    if(prod(size(Q)) < QUAD_EIG_MAX_SIZE) %#ok<PSIZE>
        if(isinf(rl))
            chkQCConvex(Q,i,'upper bound constraint',solver);
        elseif(isinf(ru))
            chkQCConvex(-Q,i,'lower bound constraint',solver);
        elseif(~isinf(rl) && ~isinf(ru)) %bit of waste really...
            chkQCConvex(Q,i,'upper bound constraint',solver);
            chkQCConvex(-Q,i,'lower bound constraint',solver);
        else
            error('Quadratic Constraint %d has infinite lower and upper bounds!');
        end
    end
end

function chkQCConvex(Q,i,str,solver)
try
   e = eig(Q);
catch %#ok<CTCH>
   e = eigs(Q);
end
if(any(e < 0)) %assuming symmetric...
    if(strcmpi(solver,'cplex'))
        error('Quadratic Constraint Q matrix in cell %d (%s) is not positive semi-definite! CPLEX only solves problems with non-convex objective terms. Please use a global solver such as SCIP to solve this problem.',i,str);
    else
        error('Quadratic Constraint Q matrix in cell %d (%s) is not positive semi-definite! Please specify a global solver such as SCIP (or general nonlinear solver such as IPOPT) to solve this problem.',i,str);
    end
end

function sdout = chkSDCone(sdcone,i,ndec, warn)
%Check size, convert to sparse([C(:) A0(:) A1(:) ... ])
[m, n] = size(sdcone);
%If n = ndec+1, already in column format (or only 1x1 matrices, same diff)
if((n-1) == ndec)
    %Check we can make a square matrix from number of rows (assume full symmetric)
    if(floor(sqrt(m)) ~= sqrt(m))
        error('Constraint Semidefinite Cone %d (assumed in column format [C(:) A0(:) ...] or [F0(:) F1(:) ...]) does not contain the correct number of rows to form a square matrix',i);
    end
    %All ok
    sdout = sparse(sdcone);
else
    %Check for full matrices [C A0 A1 ...]
    if(m*(ndec+1) ~= n)
        error('Constraint Semidefinite Cone %d (assumed in full matrix format [C A0 ...] or [F0 F1 ..]) does not contain the correct number of columns',i);
    end
    %Convert to sparse column format
    Aind = m;
    sdout = zeros(m^2,ndec+1);
    C = sdcone(:,1:m);
    sdout(:,1) = C(:);
    for i = 1:ndec
        A = sdcone(:,Aind+1:Aind+m);
        sdout(:,i+1) = A(:);
        Aind = Aind + m;
    end
    sdout = sparse(sdout);
end
%Check symmetry
m = size(sdout,1);
for j = 1:ndec+1
    C = reshape(sdout(:,j),sqrt(m),sqrt(m));
    if(abs(sum(sum(triu(C) - tril(C)'))) > 1e-6)
        %Force Symmetric (bad really)
        C = (C+C')./2;
        sdout(:,j) = C(:); %#ok<SPRIX>
        %Present Warning
        if(warn)
            if(j==1)
                optiwarn('opti:symm','Constraint Semidefinite Cone %d [matrix C (F0)] is not symmetric! Forcing to Symmetric using (C+C'')/2 but this may not be Intended.',i);
            else
                optiwarn('opti:symm','Constraint Semidefinite Cone %d [matrix A%d (F%d)] is not symmetric! Forcing to Symmetric using (A+A'')/2 but this may not be Intended.',i,j-1,j-1);
            end
        end
    end
end

function checkDerivs(prob,warn)
%Check User Supplied Derivatives

if(any(strcmpi(prob.type,{'LP','MILP','BILP','QP','MIQP','QCQP','MIQCQP','SDP','MISDP','DNLS'}))) %dnls handled internally
    return;
end

chkGrad = false;
chkJac = false;
chkHess = false;
genX0 = false;

%Read or Generate x0
if(isempty(prob.x0) || any(isnan(prob.x0)))
    if(prob.sizes.ndec == 0)
        error('OPTI cannot check your derivatives without knowing the size of the problem. Please supply x0 or ndec to the OPTI constructor.');
    end
    x0 = rand(prob.sizes.ndec,1);
    if(~isempty(prob.lb))
        ind = x0 < prob.lb;
        x0(ind) = prob.lb(ind);
    end
    if(~isempty(prob.ub))
        ind = x0 > prob.ub;
        x0(ind) = prob.ub(ind);
    end
    genX0 = true;
else
    x0 = prob.x0;
end

%Check Gradient?
if(~isempty(prob.f) && isa(prob.f,'function_handle') && ~prob.numdif.grad)
    chkGrad = true;
end
%Check Jacobian?
if(~isempty(prob.nljac) && isa(prob.nljac,'function_handle') && ~prob.numdif.jac)
    chkJac = true;
end
%Check Hessian?
if(~isempty(prob.H) && isa(prob.H,'function_handle') && ~prob.numdif.hess && chkGrad)
    chkHess = true;
end

if(genX0 && (chkGrad || chkJac || chkHess))
    if(warn>1), optiwarn('OPTI:Nox0','OPTI is randomly choosing a point a to test the user supplied derivatives. Supply a starting point as x0 to OPTI to select this point manually.'); end
end

%Do Actual Checking
if(chkGrad), optiDerCheck(prob.fun,prob.f,x0,'Objective Gradient',warn); end
if(chkJac), optiDerCheck(prob.nlcon,prob.nljac,x0,'Constraint Jacobian',warn,prob.nljacstr()); end
if(chkHess)
    if(chkJac)
        optiHessCheck(prob.H,prob.f,prob.nljac,x0,'Hessian of the Lagrangian',warn,prob.Hstr());
    else
        optiHessCheck(prob.H,prob.f,prob.nljac,x0,'Hessian',warn,prob.Hstr());
    end
end

% R2017b+ - See https://au.mathworks.com/help/optim/problem-based-lp-milp.html
function [prob,opts] = convertMatlabOptimProbToOptiProb(inProb, inOpts)

if (nargin < 2), inOpts = []; end
opts = optiset(inOpts);

try
    mlProb = prob2struct(inProb);
catch ME
    error('Error calling prob2struct on Optimization Problem - check the original problem for errors\n\n%s',ME.message);
end

% Convert known fields
prob = optiprob('f', mlProb.f, ...
                'objbias', mlProb.f0, ...
                'ineq', mlProb.Aineq, mlProb.bineq, ...
                'eq', mlProb.Aeq, mlProb.beq, ...
                'bounds', mlProb.lb, mlProb.ub, ...
                'int', mlProb.intcon);

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef MimPlugin < handle
    % MimPlugin. Base class for a Plugin used by MIM.
    %
    %     Plugins are classes you create to run your own routines from the
    %     MIM toolkit. When you create a plugin, it can automatically
    %     appear as a button in compatible guis, and is available to scripts and other
    %     plugins which use the MIM toolkit. The framework automatically
    %     handles result caching, dependency tracking and preview thumbnail
    %     generation.
    %
    %     To run a plugin from your own code, first create a suitable
    %     MimMain object, usually by using a specific application such as
    %     PTKMain for PTK. Then call CreateDatasetFromInfo() to create a MimDataset for the image files
    %     you wish to run the plugin with. Then call GetResult() on this dataset
    %     interface to run the plugin or fetch a previously cached result.
    %
    %     Example
    %     -------
    %     Replace <image path> and <filenames> with the path and filenames
    %     to your image data, and MyPluginName with the name of the plugin to
    %     run.
    %
    %         image_info = PTKImageInfo( <image path>, <filenames>, [], [], [], []);
    %         ptk = PTKMain();
    %         dataset = ptk.CreateDatasetFromInfo(image_info);
    %         results = dataset.GetResult('MyPluginName');
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    % You must set each of these properties in your plugin class
    properties (Abstract = true)
        % Should normally be set to true.
        % Specifies if the results of this plugin can be cached on disk. If you
        % set this to false, the results will never be cached and the plugin
        % will have to be run each time its result is requested. This will save
        % disk space but could substantially increase execution time.
        AllowResultsToBeCached

        % Should normally be set to false.
        % Set to true if you want to force the plugin to be run every time its
        % result is requested, ignoring any results in the disk cache. One
        % reason to set this to true is during plugin development, where you are
        % modifying the code and want to test your changes.
        AlwaysRunPlugin

        % Should normally be set to 'ReplaceOverlay'.
        % This specifies how the plugin result will be displayed to the user
        % when using the gui. Allowable values are:
        % 'ReplaceOverlay' - the result is displayed as an image overlay
        % 'ReplaceImage'   - the result replaces the image. Normally you would
        %                    only do this when creating a new image context
        %                    (i.e. region of interest)
        % 'DoNothing'      - do not display any results. This option is
        %                    appropriate if your plugin displays its results in
        %                    its own window
        PluginType

        % Set to force the plugin to appear in a particular panel name in the gui.
        % If you leave blank, the panel name will be the subdirectory under
        % Plugins where the plugin is found
        Category

        % Should normally be set to false.
        % If true, this plugin will not have a button in the gui window
        HidePluginInDisplay

        % Should normally be set to 6.
        % The width (in units defined by the gui) of the button for this plugin
        ButtonWidth

        % Should normally be set to 2.
        % The height (in units defined by the gui) of the button for this plugin
        ButtonHeight

        % The text to appear in the gui when the user hovers the mouse over the
        % button for this plugin
        ToolTip

        % The text to appear on the button in the gui. Also used in reporting of
        % progress, so the text should make sense when displayed as "Calculating
        % <button text>", e.g. 'Airways by Lobe' will display as "Calculating
        % Airways by Lobe" in the progress dialog. Note that HTML tags are
        % allowed - use <BR> for a newline in the button text
        ButtonText

        % Should normaly be set to true.
        % Specifies in a preview thumnail should be generated from the plugin
        % result after the plugin has been run. This thumbnail will be displaued
        % in the button itself. Set this to false if the plugin does not
        % generate an image.
        GeneratePreview

        % Should normally be set to false.
        % Specifies whether a MIP will be used to generate the preview
        % thumbnail. Set to true for small structures such as airways which are
        % best visualised with the "flatten" plugin.
        FlattenPreviewImage
    end

    % Optional properties
    %   PluginInterfaceVersion - specifies the version of the plugin interface
    %   that your plugin implements

    methods (Abstract, Static)

        % Generates the results for this plugin.
        % The results can be a PTKImage or
        % any kind of Matlab structure. Plugins should generally return an
        % output image as their result, but some plugins (such as PTKAirways)
        % generate information which cannot be stored in a simple image. For
        % these types of plugin, you should return a custom data structure
        % containing your results, and then override the function
        % GenerateImageFromResults below to produce an illustrative output
        % image.
        results = RunPlugin(dataset, context, reporting)

    end

    methods (Static)

        function results = GenerateImageFromResults(results, ~, ~)
            % Generates a results image from the results returned by RunPlugin.
            % If the result is itself an image, there is no need to override this
            % function as this image will be returned. If the result is a data
            % structure, this function should generate an illustrative output image
            % based on the results. For example, the PTKAirways plugin returns a
            % heirarchical data structure describing the airway tree. This functin
            % is then used to turn the structure into an image.
        end

        function result = GetEditedResult(result, edited_result, ~)
            % Modifies a plugin result based on manual editing of the output image.
            % The default behaviour is that the edited result replaces the
            % result, provided they are both of the same type.
            % However, plugins may choose to modify the original result based on
            % an edited result, in which case the plugin must override this
            % method and define its desired behaviour.
            if isempty(result) || strcmp(class(result), class(edited_result))
                result = edited_result;
            end
        end

        function result = GenerateDefaultEditedResultFollowingFailure(dataset, context, reporting)
            % Defines the behaviour if the plugin execution fails with an
            % exception. If this method returns an empty result, then the
            % exception will be rethrown and propagated down to the caller,
            % which is the default behaviour. A plugin may override this
            % method and return a default edited result, in which case the
            % user will be prompted if they wish to accept the edited
            % result.
            %
            result = [];
        end

        function enabled = IsEnabled(mim_gui_app)
            enabled = true;
        end

        function is_selected = IsSelected(mim_gui_app)
            is_selected = false;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function A = randomize_graph_partial_und(A,B,maxswap)
% RANDOMIZE_GRAPH_PARTIAL_UND    Swap edges with preserved degree sequence
%
%   A = RANDOMIZE_GRAPH_PARTIAL_UND(A,B,MAXSWAP) takes adjacency matrices A
%   and B and attempts to randomize matrix A by performing MAXSWAP
%   rewirings. The rewirings will avoid any spots where matrix B is
%   nonzero.
%
%   Inputs:       A,      undirected adjacency matrix
%                 B,      edges to avoid
%           MAXSWAP,      number of rewirings
%
%   Outputs:      A,      randomized matrix
%
%   Richard Betzel, Indiana University, 2013
%
%   Notes:
%   1. Based on the script randmio_und.m.
%   2. Graph may become disconnected as a result of rewiring. Always
%      important to check.
%   3. A can be weighted, though the weighted degree sequence will not be
%      preserved.
%

[i,j] = find(triu(A,1));
m = length(i);
nswap = 0;
while nswap < maxswap
    while 1
        e1 = randi(m); e2 = randi(m);
        while e2 == e1
            e2 = randi(m);
        end
        a = i(e1); b = j(e1);
        c = i(e2); d = j(e2);
        if all(a~=[c,d]) && all(b~=[c,d])
            break
        end
    end
    if rand > 0.5
        i(e2) = d; j(e2) = c;
        c = i(e2); d = j(e2);
    end
    if ~(A(a,d) || A(c,b) || B(a,d) || B(c,b))
        A(a,d) = A(a,b); A(a,b) = 0;
        A(d,a) = A(b,a); A(b,a) = 0;
        A(c,b) = A(c,d); A(c,d) = 0;
        A(b,c) = A(d,c); A(d,c) = 0;
        j(e1) = d;
        j(e2) = b;
        nswap = nswap + 1;
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function h=error_ellipse(varargin)
% ERROR_ELLIPSE - plot an error ellipse, or ellipsoid, defining confidence region
%    ERROR_ELLIPSE(C22) - Given a 2x2 covariance matrix, plot the
%    associated error ellipse, at the origin. It returns a graphics handle
%    of the ellipse that was drawn.
%
%    ERROR_ELLIPSE(C33) - Given a 3x3 covariance matrix, plot the
%    associated error ellipsoid, at the origin, as well as its projections
%    onto the three axes. Returns a vector of 4 graphics handles, for the
%    three ellipses (in the X-Y, Y-Z, and Z-X planes, respectively) and for
%    the ellipsoid.
%
%    ERROR_ELLIPSE(C,MU) - Plot the ellipse, or ellipsoid, centered at MU,
%    a vector whose length should match that of C (which is 2x2 or 3x3).
%
%    ERROR_ELLIPSE(...,'Property1',Value1,'Name2',Value2,...) sets the
%    values of specified properties, including:
%      'C' - Alternate method of specifying the covariance matrix
%      'mu' - Alternate method of specifying the ellipse (-oid) center
%      'conf' - A value betwen 0 and 1 specifying the confidence interval.
%        the default is 0.5 which is the 50% error ellipse.
%      'scale' - Allow the plot the be scaled to difference units.
%      'style' - A plotting style used to format ellipses.
%      'clip' - specifies a clipping radius. Portions of the ellipse, -oid,
%        outside the radius will not be shown.
%
%    NOTES: C must be positive definite for this function to work properly.

default_properties = struct(...
  'C', [], ... % The covaraince matrix (required)
  'mu', [], ... % Center of ellipse (optional)
  'conf', 0.5, ... % Percent confidence/100
  'scale', 1, ... % Scale factor, e.g. 1e-3 to plot m as km
  'style', '', ...  % Plot style
  'clip', inf); % Clipping radius

if length(varargin) >= 1 & isnumeric(varargin{1})
  default_properties.C = varargin{1};
  varargin(1) = [];
end

if length(varargin) >= 1 & isnumeric(varargin{1})
  default_properties.mu = varargin{1};
  varargin(1) = [];
end

if length(varargin) >= 1 & isnumeric(varargin{1})
  default_properties.conf = varargin{1};
  varargin(1) = [];
end

if length(varargin) >= 1 & isnumeric(varargin{1})
  default_properties.scale = varargin{1};
  varargin(1) = [];
end

if length(varargin) >= 1 & ~ischar(varargin{1})
  error('Invalid parameter/value pair arguments.')
end

prop = getopt(default_properties, varargin{:});
C = prop.C;

if isempty(prop.mu)
  mu = zeros(length(C),1);
else
  mu = prop.mu;
end

conf = prop.conf;
scale = prop.scale;
style = prop.style;

if conf <= 0 | conf >= 1
  error('conf parameter must be in range 0 to 1, exclusive')
end

[r,c] = size(C);
if r ~= c | (r ~= 2 & r ~= 3)
  error(['Don''t know what to do with ',num2str(r),'x',num2str(c),' matrix'])
end

x0=mu(1);
y0=mu(2);

% Compute quantile for the desired percentile
k = sqrt(qchisq(conf,r)); % r is the number of dimensions (degrees of freedom)

hold_state = get(gca,'nextplot');

if r==3 & c==3
  z0=mu(3);

  % Make the matrix has positive eigenvalues - else it's not a valid covariance matrix!
  if any(eig(C) <=0)
    error('The covariance matrix must be positive definite (it has non-positive eigenvalues)')
  end

  % C is 3x3; extract the 2x2 matricies, and plot the associated error
  % ellipses. They are drawn in space, around the ellipsoid; it may be
  % preferable to draw them on the axes.
  Cxy = C(1:2,1:2);
  Cyz = C(2:3,2:3);
  Czx = C([3 1],[3 1]);

  [x,y,z] = getpoints(Cxy,prop.clip);
  h1=plot3(x0+k*x,y0+k*y,z0+k*z,prop.style);hold on
  [y,z,x] = getpoints(Cyz,prop.clip);
  h2=plot3(x0+k*x,y0+k*y,z0+k*z,prop.style);hold on
  [z,x,y] = getpoints(Czx,prop.clip);
  h3=plot3(x0+k*x,y0+k*y,z0+k*z,prop.style);hold on

  [eigvec,eigval] = eig(C);

  [X,Y,Z] = ellipsoid(0,0,0,1,1,1);
  XYZ = [X(:),Y(:),Z(:)]*sqrt(eigval)*eigvec';

  X(:) = scale*(k*XYZ(:,1)+x0);
  Y(:) = scale*(k*XYZ(:,2)+y0);
  Z(:) = scale*(k*XYZ(:,3)+z0);
  h4=surf(X,Y,Z);
  colormap gray
  alpha(0.3)
  camlight
  if nargout
    h=[h1 h2 h3 h4];
  end
elseif r==2 & c==2
  % Make the matrix has positive eigenvalues - else it's not a valid covariance matrix!
  if any(eig(C) <=0)
    error('The covariance matrix must be positive definite (it has non-positive eigenvalues)')
  end

  [x,y,z] = getpoints(C,prop.clip);
  h1=plot(scale*(x0+k*x),scale*(y0+k*y),prop.style);
  set(h1,'zdata',z+1)
  if nargout
    h=h1;
  end
else
  error('C (covaraince matrix) must be specified as a 2x2 or 3x3 matrix)')
end
%axis equal

set(gca,'nextplot',hold_state);

%---------------------------------------------------------------
% getpoints - Generate x and y points that define an ellipse, given a 2x2
%   covariance matrix, C. z, if requested, is all zeros with same shape as
%   x and y.
function [x,y,z] = getpoints(C,clipping_radius)

n=100; % Number of points around ellipse
p=0:pi/n:2*pi; % angles around a circle

[eigvec,eigval] = eig(C); % Compute eigen-stuff
xy = [cos(p'),sin(p')] * sqrt(eigval) * eigvec'; % Transformation
x = xy(:,1);
y = xy(:,2);
z = zeros(size(x));

% Clip data to a bounding radius
if nargin >= 2
  r = sqrt(sum(xy.^2,2)); % Euclidian distance (distance from center)
  x(r > clipping_radius) = nan;
  y(r > clipping_radius) = nan;
  z(r > clipping_radius) = nan;
end

%---------------------------------------------------------------
function x=qchisq(P,n)
% QCHISQ(P,N) - quantile of the chi-square distribution.
if nargin<2
  n=1;
end

s0 = P==0;
s1 = P==1;
s = P>0 & P<1;
x = 0.5*ones(size(P));
x(s0) = -inf;
x(s1) = inf;
x(~(s0|s1|s))=nan;

for ii=1:14
  dx = -(pchisq(x(s),n)-P(s))./dchisq(x(s),n);
  x(s) = x(s)+dx;
  if all(abs(dx) < 1e-6)
    break;
  end
end

%---------------------------------------------------------------
function F=pchisq(x,n)
% PCHISQ(X,N) - Probability function of the chi-square distribution.
if nargin<2
  n=1;
end
F=zeros(size(x));

if rem(n,2) == 0
  s = x>0;
  k = 0;
  for jj = 0:n/2-1;
    k = k + (x(s)/2).^jj/factorial(jj);
  end
  F(s) = 1-exp(-x(s)/2).*k;
else
  for ii=1:numel(x)
    if x(ii) > 0
      F(ii) = quadl(@dchisq,0,x(ii),1e-6,0,n);
    else
      F(ii) = 0;
    end
  end
end

%---------------------------------------------------------------
function f=dchisq(x,n)
% DCHISQ(X,N) - Density function of the chi-square distribution.
if nargin<2
  n=1;
end
f=zeros(size(x));
s = x>=0;
f(s) = x(s).^(n/2-1).*exp(-x(s)/2)./(2^(n/2)*gamma(n/2));

%---------------------------------------------------------------
function properties = getopt(properties,varargin)
%GETOPT - Process paired optional arguments as 'prop1',val1,'prop2',val2,...
%
%   getopt(properties,varargin) returns a modified properties structure,
%   given an initial properties structure, and a list of paired arguments.
%   Each argumnet pair should be of the form property_name,val where
%   property_name is the name of one of the field in properties, and val is
%   the value to be assigned to that structure field.
%
%   No validation of the values is performed.
%
% EXAMPLE:
%   properties = struct('zoom',1.0,'aspect',1.0,'gamma',1.0,'file',[],'bg',[]);
%   properties = getopt(properties,'aspect',0.76,'file','mydata.dat')
% would return:
%   properties =
%         zoom: 1
%       aspect: 0.7600
%        gamma: 1
%         file: 'mydata.dat'
%           bg: []
%
% Typical usage in a function:
%   properties = getopt(properties,varargin{:})

% Process the properties (optional input arguments)
prop_names = fieldnames(properties);
TargetField = [];
for ii=1:length(varargin)
  arg = varargin{ii};
  if isempty(TargetField)
    if ~ischar(arg)
      error('Propery names must be character strings');
    end
    f = find(strcmp(prop_names, arg));
    if length(f) == 0
      error('%s ',['invalid property ''',arg,'''; must be one of:'],prop_names{:});
    end
    TargetField = arg;
  else
    % properties.(TargetField) = arg; % Ver 6.5 and later only
    properties = setfield(properties, TargetField, arg); % Ver 6.1 friendly
    TargetField = '';
  end
end
if ~isempty(TargetField)
  error('Property names and values must be specified in pairs.');
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function planC = flipToMatchPACS(scanNumV,dim,planC)
% flipToMatchPACS(planC,scanNum,dim)
% ---- INPUTS ----
% scanNumV: Indices of scans to be flipped
% dim     : May be 'lr', 'ap', 'si'
% planC   : CERR archive
% ---------------
% AI 9/20/17

indexS = planC{end};

for i = 1:length(scanNumV)
    nSlice = length([planC{indexS.scan}(scanNumV(i)).scanInfo]);
    switch lower(dim)

        case 'lr'
            if ~isfield(planC{indexS.scan}(scanNumV(i)).scanInfo(1),'LRflippedToMatchPACS')
                [planC{indexS.scan}(scanNumV(i)).scanInfo(:).LRflippedToMatchPACS] = deal('');
            end
            planC = flipAlongX(scanNumV(i),planC);
            currentFlagV = [planC{indexS.scan}(scanNumV(i)).scanInfo(:).LRflippedToMatchPACS];
            if isempty(currentFlagV)
                currentFlagV = false(nSlice,1);
            end
            flipFlagC = num2cell(~currentFlagV);
            [planC{indexS.scan}(scanNumV(i)).scanInfo(:).LRflippedToMatchPACS] = flipFlagC{:};

        case 'ap'
            if ~isfield(planC{indexS.scan}(scanNumV(i)).scanInfo(1),'APflippedToMatchPACS')
                [planC{indexS.scan}(scanNumV(i)).scanInfo(:).APflippedToMatchPACS] = deal('');
            end
            planC = flipAlongY(scanNumV(i),planC);
            currentFlagV = [planC{indexS.scan}(scanNumV(i)).scanInfo(:).APflippedToMatchPACS];
            if isempty(currentFlagV)
                currentFlagV = false(nSlice,1);
            end
            [planC{indexS.scan}(scanNumV(i)).scanInfo(:).APflippedToMatchPACS] = ~currentFlagV;

        case 'si'
            if ~isfield(planC{indexS.scan}(scanNumV(i)).scanInfo(1),'SIflippedToMatchPACS')
                [planC{indexS.scan}(scanNumV(i)).scanInfo(:).SIflippedToMatchPACS] = deal('');
            end
            planC = flipAlongZ(scanNumV(i),planC);
            currentFlagV = [planC{indexS.scan}(scanNumV(i)).scanInfo(:).SIflippedToMatchPACS];
            if isempty(currentFlagV)
                currentFlagV = false(nSlice,1);
            end
            [planC{indexS.scan}(scanNumV(i)).scanInfo(:).SIflippedToMatchPACS] = ~currentFlagV;

        otherwise
            CERRString('Invalid input. ''dim'' must be ''lr'' ,''ap'' or ''si'' ','console');
            return
    end

end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% Maximum volume inscribed ellipsoid in a polyhedron
% Section 8.4.1, Boyd & Vandenberghe "Convex Optimization"
% Original version by Lieven Vandenberghe
% Updated for CVX by Almir Mutapcic - Jan 2006
% (a figure is generated)
%
% We find the ellipsoid E of maximum volume that lies inside of
% a polyhedra C described by a set of linear inequalities.
%
% C = { x | a_i^T x <= b_i, i = 1,...,m } (polyhedra)
% E = { Bu + d | || u || <= 1 } (ellipsoid)
%
% This problem can be formulated as a log det maximization
% which can then be computed using the det_rootn function, ie,
%     maximize     log det B
%     subject to   || B a_i || + a_i^T d <= b,  for i = 1,...,m

% problem data
n = 2;
px = [0 .5 2 3 1];
py = [0 1 1.5 .5 -.5];
m = size(px,2);
pxint = sum(px)/m; pyint = sum(py)/m;
px = [px px(1)];
py = [py py(1)];

% generate A,b
A = zeros(m,n); b = zeros(m,1);
for i=1:m
  A(i,:) = null([px(i+1)-px(i) py(i+1)-py(i)])';
  b(i) = A(i,:)*.5*[px(i+1)+px(i); py(i+1)+py(i)];
  if A(i,:)*[pxint; pyint]-b(i)>0
    A(i,:) = -A(i,:);
    b(i) = -b(i);
  end
end

% formulate and solve the problem
cvx_begin
    variable B(n,n) symmetric
    variable d(n)
    maximize( det_rootn( B ) )
    subject to
       for i = 1:m
           norm( B*A(i,:)', 2 ) + A(i,:)*d <= b(i); %#ok
       end
cvx_end

% make the plots
noangles = 200;
angles   = linspace( 0, 2 * pi, noangles );
ellipse_inner  = B * [ cos(angles) ; sin(angles) ] + d * ones( 1, noangles );
ellipse_outer  = 2*B * [ cos(angles) ; sin(angles) ] + d * ones( 1, noangles );

clf
plot(px,py)
hold on
plot( ellipse_inner(1,:), ellipse_inner(2,:), 'r--' );
plot( ellipse_outer(1,:), ellipse_outer(2,:), 'r--' );
axis square
axis off
hold off

author: yu-jiang

repo: radpbook

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


text:
function [variables, index] = findGUBGroup(p,s)
variables = [];index = [];
for i = find((p.knapsack.type == 3) | (p.knapsack.type == 5))
    if all(ismember(s,p.knapsack.variables{i}))
        variables = p.knapsack.variables{i};
        index = i;
        return
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [v,n]=cutvar(t,j)
%CUTVAR Variables used for branches in decision tree.
%   V=CUTVAR(T) returns an N-element cell array V of the names of the
%   variables used for branching in each node of the tree T, where N
%   is the number of nodes in the tree.  These variables are sometimes
%   known as cut variables.  For non-branch (leaf) nodes, V contains
%   an empty string.
%
%   V=CUTVAR(T,J) takes an array J of node numbers and returns the cut
%   variables for the specified nodes.
%
%   [V,NUM]=CUTVAR(...) also returns the number NUM of each variable.
%
%   See also CLASSREGTREE, CLASSREGTREE/NUMNODES, CLASSREGTREE/CHILDREN.

%   Copyright 2006-2007 The MathWorks, Inc.
%   $Revision: 1.1.6.2 $  $Date: 2007/02/15 21:48:07 $

if nargin>=2 && ~validatenodes(t,j)
    error('stats:classregtree:parent:InvalidNode',...
          'J must be an array of node numbers or a logical array of the proper size.');
end

% Get variable number, taking care that negative values may appear to
% indicate categorical cuts
if nargin<2
    n = abs(t.var);
else
    n = abs(t.var(j,:));
end

% Get variable name
v = repmat({''},numel(n),1);
mask = n>0;
v(mask) = t.names(n(mask));

author: Cloud-CV

repo: object-proposals


text:
% this is an experimentaal high-performance EPI sequence
% which uses split gradients to overlap blips with the readout
% gradients combined with ramp-samping

seq=mr.Sequence();         % Create a new sequence object
fov=250e-3; Nx=64; Ny=64;  % Define FOV and resolution
thickness=3e-3;            % slice thinckness
Nslices=3;
TE=40e-3;

pe_enable=1;               % a flag to quickly disable phase encoding (1/0) as needed for the delay calibration
ro_os=1;                   % oversampling factor (in contrast to the product sequence we don't really need it)
readoutTime=4.2e-4;        % this controls the readout bandwidth
partFourierFactor=0.75;    % partial Fourier factor: 1: full sampling 0: start with ky=0

tRFex=2e-3;
tRFref=2e-3;
spoilFactor=1.5;             % spoiling gradient around the pi-pulse

% Set system limits
lims = mr.opts('MaxGrad',32,'GradUnit','mT/m',...
    'MaxSlew',130,'SlewUnit','T/m/s',...
    'rfRingdownTime', 30e-6, 'rfDeadtime', 100e-6);

% Create fat-sat pulse
B0=2.89; % 1.5 2.89 3.0
sat_ppm=-3.45;
sat_freq=sat_ppm*1e-6*B0*lims.gamma;
rf_fs = mr.makeGaussPulse(110*pi/180,'system',lims,'Duration',8e-3,...
    'bandwidth',abs(sat_freq),'freqOffset',sat_freq);
rf_fs.phaseOffset=-2*pi*rf_fs.freqOffset*mr.calcRfCenter(rf_fs); % compensate for the frequency-offset induced phase
gz_fs = mr.makeTrapezoid('z',lims,'delay',mr.calcDuration(rf_fs),'Area',1/1e-4); % spoil up to 0.1mm

% Create 90 degree slice selection pulse and gradient
[rf, gz, gzReph] = mr.makeSincPulse(pi/2,'system',lims,'Duration',tRFex,...
    'SliceThickness',thickness,'apodization',0.5,'timeBwProduct',4);

% Create 90 degree slice refocusing pulse and gradients
[rf180, gz180] = mr.makeSincPulse(pi,'system',lims,'Duration',tRFref,...
    'SliceThickness',thickness,'apodization',0.5,'timeBwProduct',4,'PhaseOffset',pi/2,'use','refocusing');
% [~, gzr_t, gzr_a]=mr.makeExtendedTrapezoidArea('z',gz180.amplitude,0,-gzReph.area+0.5*gz180.amplitude*gz180.fallTime,lims);
% gz180n=mr.makeExtendedTrapezoid('z','system',lims,'times',[0 gz180.riseTime gz180.riseTime+gz180.flatTime+gzr_t]+gz180.delay, 'amplitudes', [0 gz180.amplitude gzr_a]);
[~, gzr1_t, gzr1_a]=mr.makeExtendedTrapezoidArea('z',0,gz180.amplitude,spoilFactor*gz.area,lims);
[~, gzr2_t, gzr2_a]=mr.makeExtendedTrapezoidArea('z',gz180.amplitude,0,-gzReph.area+spoilFactor*gz.area,lims);
if gz180.delay>(gzr1_t(4)-gz180.riseTime)
    gz180.delay=gz180.delay-(gzr1_t(4)-gz180.riseTime);
else
    rf180.delay=rf180.delay+(gzr1_t(4)-gz180.riseTime)-gz180.delay;
    gz180.delay=0;
end
gz180n=mr.makeExtendedTrapezoid('z','system',lims,'times',[gzr1_t gzr1_t(4)+gz180.flatTime+gzr2_t]+gz180.delay, 'amplitudes', [gzr1_a gzr2_a]);

% define the output trigger to play out with every slice excitatuion
trig=mr.makeDigitalOutputPulse('osc0','duration', 100e-6); % possible channels: 'osc0','osc1','ext1'

% Define other gradients and ADC events
deltak=1/fov;
kWidth = Nx*deltak;

% Phase blip in shortest possible time
blip_dur = ceil(2*sqrt(deltak/lims.maxSlew)/10e-6/2)*10e-6*2; % we round-up the duration to 2x the gradient raster time
% the split code below fails if this really makes a trpezoid instead of a triangle...
gy = mr.makeTrapezoid('y',lims,'Area',-deltak,'Duration',blip_dur); % we use negative blips to save one k-space line on our way towards the k-space center
%gy = mr.makeTrapezoid('y',lims,'amplitude',deltak/blip_dur*2,'riseTime',blip_dur/2, 'flatTime', 0);

% readout gradient is a truncated trapezoid with dead times at the beginnig
% and at the end each equal to a half of blip_dur
% the area between the blips should be defined by kWidth
% we do a two-step calculation: we first increase the area assuming maximum
% slewrate and then scale down the amlitude to fix the area
extra_area=blip_dur/2*blip_dur/2*lims.maxSlew; % check unit!;
gx = mr.makeTrapezoid('x',lims,'Area',kWidth+extra_area,'duration',readoutTime+blip_dur);
actual_area=gx.area-gx.amplitude/gx.riseTime*blip_dur/2*blip_dur/2/2-gx.amplitude/gx.fallTime*blip_dur/2*blip_dur/2/2;
gx.amplitude=gx.amplitude/actual_area*kWidth;
gx.area = gx.amplitude*(gx.flatTime + gx.riseTime/2 + gx.fallTime/2);
gx.flatArea = gx.amplitude*gx.flatTime;

% calculate ADC
% we use ramp sampling, so we have to calculate the dwell time and the
% number of samples, which are will be qite different from Nx and
% readoutTime/Nx, respectively.
adcDwellNyquist=deltak/gx.amplitude/ro_os;
% round-down dwell time to 100 ns
adcDwell=floor(adcDwellNyquist*1e7)*1e-7;
adcSamples=floor(readoutTime/adcDwell/4)*4; % on Siemens the number of ADC samples need to be divisible by 4
% MZ: no idea, whether ceil,round or floor is better for the adcSamples...
adc = mr.makeAdc(adcSamples,'Dwell',adcDwell,'Delay',blip_dur/2);
% realign the ADC with respect to the gradient
time_to_center=adc.dwell*((adcSamples-1)/2+0.5); % I've been told that Siemens samples in the center of the dwell period
adc.delay=round((gx.riseTime+gx.flatTime/2-time_to_center)*1e6)*1e-6; % we adjust the delay to align the trajectory with the gradient. We have to aligh the delay to 1us
% this rounding actually makes the sampling points on odd and even readouts
% to appear misalligned. However, on the real hardware this misalignment is
% much stronger anyways due to the grdient delays

% FOV positioning requires alignment to grad. raster... -> TODO

% split the blip into two halves and produnce a combined synthetic gradient
gy_parts = mr.splitGradientAt(gy, blip_dur/2, lims);
[gy_blipup, gy_blipdown, ~]=mr.align('right',gy_parts(1),'left',gy_parts(2),gx);
gy_blipdownup=mr.addGradients({gy_blipdown, gy_blipup}, lims);

% pe_enable support
gy_blipup.waveform=gy_blipup.waveform*pe_enable;
gy_blipdown.waveform=gy_blipdown.waveform*pe_enable;
gy_blipdownup.waveform=gy_blipdownup.waveform*pe_enable;

% phase encoding and partial Fourier

Ny_pre=round(partFourierFactor*Ny/2-1); % PE steps prior to ky=0, excluding the central line
Ny_post=round(Ny/2+1); % PE lines after the k-space center including the central line
Ny_meas=Ny_pre+Ny_post;

% Pre-phasing gradients
gxPre = mr.makeTrapezoid('x',lims,'Area',-gx.area/2);
gyPre = mr.makeTrapezoid('y',lims,'Area',Ny_pre*deltak);
[gxPre,gyPre]=mr.align('right',gxPre,'left',gyPre);
% relax the PE prepahser to reduce stimulation
gyPre = mr.makeTrapezoid('y',lims,'Area',gyPre.area,'Duration',mr.calcDuration(gxPre,gyPre));
gyPre.amplitude=gyPre.amplitude*pe_enable;

% Calculate delay times
durationToCenter = (Ny_pre+0.5)*mr.calcDuration(gx);
rfCenterInclDelay=rf.delay + mr.calcRfCenter(rf);
rf180centerInclDelay=rf180.delay + mr.calcRfCenter(rf180);
delayTE1=ceil((TE/2 - mr.calcDuration(rf,gz) + rfCenterInclDelay - rf180centerInclDelay)/lims.gradRasterTime)*lims.gradRasterTime;
delayTE2=ceil((TE/2 - mr.calcDuration(rf180,gz180n) + rf180centerInclDelay - durationToCenter)/lims.gradRasterTime)*lims.gradRasterTime;
assert(delayTE1>=0);
%assert(delayTE2>=0);
% now we merge slice refocusing, TE delay and pre-phasers into a single
% block
delayTE2=delayTE2+mr.calcDuration(rf180,gz180n);
gxPre.delay=0;
gxPre.delay=delayTE2-mr.calcDuration(gxPre);
assert(gxPre.delay>=mr.calcDuration(rf180)); % gxPre may not overlap with the RF
gyPre.delay=mr.calcDuration(rf180);
assert(mr.calcDuration(gyPre)<=mr.calcDuration(gxPre)); % gyPre may not shift the timing

% Define sequence blocks

%seq.addBlock(mr.makeDelay(1)); % older scanners like Trio may need this
                                % dummy delay to keep up with timing

for s=1:Nslices
    seq.addBlock(rf_fs,gz_fs);
    rf.freqOffset=gz.amplitude*thickness*(s-1-(Nslices-1)/2);
    rf.phaseOffset=-2*pi*rf.freqOffset*mr.calcRfCenter(rf); % compensate for the slice-offset induced phase
    rf180.freqOffset=gz180.amplitude*thickness*(s-1-(Nslices-1)/2);
    rf180.phaseOffset=pi/2-2*pi*rf180.freqOffset*mr.calcRfCenter(rf180); % compensate for the slice-offset induced phase
    seq.addBlock(rf,gz,trig);
    seq.addBlock(mr.makeDelay(delayTE1));
    seq.addBlock(rf180,gz180n,mr.makeDelay(delayTE2),gxPre,gyPre);
    for i=1:Ny_meas
        if i==1
            seq.addBlock(gx,gy_blipup,adc); % Read the first line of k-space with a single half-blip at the end
        elseif i==Ny_meas
            seq.addBlock(gx,gy_blipdown,adc); % Read the last line of k-space with a single half-blip at the beginning
        else
            seq.addBlock(gx,gy_blipdownup,adc); % Read an intermediate line of k-space with a half-blip at the beginning and a half-blip at the end
        end
        gx.amplitude = -gx.amplitude;   % Reverse polarity of read gradient
    end
end

%% check whether the timing of the sequence is correct
[ok, error_report]=seq.checkTiming;

if (ok)
    fprintf('Timing check passed successfully\n');
else
    fprintf('Timing check failed! Error listing follows:\n');
    fprintf([error_report{:}]);
    fprintf('\n');
end

%% do some visualizations

seq.plot();             % Plot sequence waveforms

% trajectory calculation
[ktraj_adc1, t_adc1, ktraj1, t_ktraj1, t_excitation1, t_refocusing1] = seq.calculateKspacePP();

% plot k-spaces
figure; plot(t_ktraj1, ktraj1'); % plot the entire k-space trajectory
hold on; plot(t_adc1,ktraj_adc1(1,:),'.'); % and sampling points on the kx-axis
figure; plot(ktraj1(1,:),ktraj1(2,:),'b'); % a 2D plot
axis('equal'); % enforce aspect ratio for the correct trajectory display
hold;plot(ktraj_adc1(1,:),ktraj_adc1(2,:),'r.'); % plot the sampling points

%% prepare the sequence output for the scanner
seq.setDefinition('FOV', [fov fov thickness]);
seq.setDefinition('Name', 'epi');

seq.write('epise_rs.seq');

% seq.install('siemens');

% seq.sound(); % simulate the seq's tone

%% very optional slow step, but useful for testing during development e.g. for the real TE, TR or for staying within slewrate limits

rep = seq.testReport;
fprintf([rep{:}]); % as for January 2019 TR calculation fails for fat-sat

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq