text:
###############################################
"""
$(TYPEDEF)

Structure containing the parameters of the 2D Beta functions.

# Fields

$(TYPEDFIELDS)

"""
struct ScaledBeta2DParams
    "number of model parameters"
    nummodpar::Integer
    "type of function on y (protein concentration) for the mode"
    modefuny::String
    "type of function on y (protein concentration) for the confidence parameter"
    konfuny::String
    "type of function on y (protein concentration) for the amplitude"
    ampfuny::String
    "lower bound for Beta function (along x, i.e. SDS concentration)"
    a::Real
    "upper bound for Beta function (along y, i.e. protein concentration)"
    b::Real
    "user defined minimum y (protein concentration)"
    ymin::Real
    "user defined maximum y (protein concentration)"
    ymax::Real
    "starting index in the vector of model parameters for mode"
    idxmode::Integer
    "starting index in the vector of model parameters for confidence parameter"
    idxkon::Integer
    "starting index in the vector of model parameters for amplitude parameter"
    idxamp::Integer

    ##----------------------------------------------------------------------
    function ScaledBeta2DParams(; modefuny,konfuny,ampfuny,a,b,ymin,ymax)
        nummodpar::Integer=0
        if modefuny=="linear"
            nummodpar+=2
            idxmode=1
        elseif modefuny=="quadratic"
            nummodpar+=3
            idxmode=1
        else
            error("ScaledBeta2DParams: wrong `modefuny`")
        end
        if konfuny=="constant"
            idxkon=nummodpar+1
            nummodpar+=1
        elseif konfuny=="linear"
            idxkon=nummodpar+1
            nummodpar+=2
        elseif konfuny=="quadratic"
            idxkon=nummodpar+1
            nummodpar+=3
        else
            error("ScaledBeta2DParams: wrong `konfuny`")
        end
        if ampfuny=="linear"
            idxamp=nummodpar+1
            nummodpar+=2
        elseif ampfuny=="quadratic"
            idxamp=nummodpar+1
            nummodpar+=3
        else
            error("ScaledBeta2DParams: wrong `ampfuny`")
        end
        return new(nummodpar,modefuny,konfuny,ampfuny,a,b,ymin,ymax,idxmode,idxkon,idxamp)
    end
    ##----------------------------------------------------------------------
end

###########################################################

"""
$(TYPEDEF)

Structure containing the observed (measured) data, i.e. the enthalpy, the protein and SDS concentration.

# Fields

$(TYPEDFIELDS)

"""
Base.@kwdef struct ITCObsData
    "Enthalpy values"
    enthalpy::Vector{<:Real}
    "Indices pointing to single experiments data"
    idxdata::Vector{UnitRange{Int64}}
    "Concentration of SDS (first column) and protein (second column)"
    sdsprotcon::Array{<:Real,2}
    "Name of the protein"
    protein::String
end

###########################################################

"""
$(TYPEDEF)

Structure containing the observed (measured) data, i.e. the enthalpy, the protein and SDS concentration.

# Fields

$(TYPEDFIELDS)

"""
struct BetaMix2D
    "structure containing the parameters of the 2D Beta functions"
    betpar::ScaledBeta2DParams
    "2D array where each column is a set of mode, confidence and amplitude parameters"
    modkonamp::Matrix{<:Real}
    "name of the protein"
    protein::String
end

##############################################################

# """
# $(TYPEDSIGNATURES)

# Set the constraints for the Newton optimization.
# """
# function setconstraints(betpar::ScaledBeta2DParams,mstart::Matrix{<:Real},
#                         lowc::Vector{<:Real},upc::Vector{<:Real})

#     lenm = size(mstart,1)
#     ncomp = size(mstart,2) #div(lenm,betpar.nummodpar)

#     lowconstr = zeros(betpar.nummodpar,ncomp)
#     upconstr = zeros(betpar.nummodpar,ncomp)

#     #@show size(lowc),size(upc)
#     for c=1:ncomp
#         #c1 = (c-1)*betpar.nummodpar
#         for i=1:betpar.nummodpar
#             lowconstr[i,c] = lowc[i]
#             upconstr[i,c] = upc[i]
#         end
#     end

#     @show lowconstr
#     @show upconstr
#     @show mstart.<=lowconstr
#     @show mstart.>=upconstr

#     ## check mstart...
#     @assert all(mstart.>=lowconstr)
#     @assert all(mstart.<=upconstr)

#     return vec(lowconstr),vec(upconstr)
# end

######################################################################

"""
$(TYPEDSIGNATURES)

Calculate the sum of the misfit for observed (measured) and calculated data and
 the misfit of the integral of enthalpy (which should be zero) at zero protein
 concentration.
"""
function misfitfunctional(betpar::ScaledBeta2DParams,dobs::ITCObsData,
                          invCd::Matrix{<:Real},mcur::Matrix{<:Real} )
                          #lagmularea::Union{Real,Nothing}=nothing )

    ## misfit to measured data for Beta functions
    misflkl = misfitbeta2D(betpar,dobs,invCd,mcur)

    # # total misfit
    # if lagmularea==nothing
    #     ## no area
         misf = misflkl
    # else
    #     ## enthalpy area
    #     misfar = misfareaenth(betpar,mcur)
    #     misf = misflkl + lagmularea*misfar
    # end

    return misf
end

###############################################################

"""
$(TYPEDSIGNATURES)

Calculate misfit between observed/measured and calculated data.
"""
function misfitbeta2D(betpar::ScaledBeta2DParams,dobs::ITCObsData,
                      invCd::Matrix{<:Real},mcur::Matrix{<:Real})

    # compute synthetic data
    # xy = dobs.sdsprotcon
    dcalc = forwmod2D(betpar,dobs.sdsprotcon,mcur)

    # calculate likelihood value
    difcalobs = dcalc.-dobs.enthalpy
    tmp1 = invCd * difcalobs
    misflkl = 0.5 * dot(difcalobs,tmp1)

  return misflkl
end

###########################################################

"""
$(TYPEDSIGNATURES)

Misfit functional for area of enthalpy (unscaled in this case): zero
 area at protein concentration equal to zero.
"""
function misfareaenth(betpar::ScaledBeta2DParams,mcur2d::Matrix{<:Real})

    ## set the constraint of zero area at protcon=0.0
    protcon = 0.0
    ncomp = size(mcur2d,2)
    integr = zeros(Real,ncomp)
    err = zeros(ncomp)

    abeta = betpar.a
    bbeta = betpar.b

    ## Using Gaussian quadrature from QuadGK instead of Cubature to make ForwardDiff work
    Npts = 500
    qtpts,weights = gauss(Npts,abeta,bbeta)

    for c=1:ncomp
        mcur = mcur2d[:,c]
        @assert length(mcur)==betpar.nummodpar

        ##---------------------------------------------------------
        # (val,err) = hquadrature(f::Function, xmin::Real, xmax::Real;
        #                     reltol=1e-8, abstol=0, maxevals=0)

        ## get the values of model parameters at y=ycur
        mode,kon,amp = getmodparbeta(betpar,mcur,protcon)

        @assert all(abeta.<=qtpts.<=bbeta)
        # if !(abeta<=mode<=bbeta)
        #     @show mcur2d
        #     betami = BetaMix2D(betpar,mcur2d,"boh")
        #     plotparamlines(betami)
        # end
        @assert abeta<=mode<=bbeta
        @assert kon>2.0

        ## gaussian quadrature
        nodes = [scaledbeta(mode,kon,abeta,bbeta,amp,x) for x in qtpts]
        integr[c] = dot( weights, nodes )

    end

    areanorm = sum(abs.(integr))  ## <<<<<===========  CHECK CHECK CHECK!!
    return areanorm
end

###########################################################

"""
$(TYPEDSIGNATURES)

Compute the forward response for given input parameters, i.e. the enthalpy
 values (2D) for given Beta functions (2D) parameters.
"""
function forwmod2D(betpar::ScaledBeta2DParams,xy::Array{<:Real,2},mcur::Matrix{<:Real})

    npts = size(xy,1)
    nummodpar = size(mcur,1)
    ncomp = size(mcur,2)

    sscbeta = zeros(Real,npts) #zeros(eltype(mcur),npts)
    # sscbeta .= zero(eltype(mcur))   #0.0 ## make sure it's zeroed

    for n=1:ncomp
        # sum all the components
        sscbeta .+= singlescaledbeta2D(betpar,xy,mcur[:,n])
    end

    return sscbeta
end

###########################################################

"""
 $(TYPEDSIGNATURES)

  2D-fied scaled Beta function for given x (SDS concentration) and y (protein concentration) values.
"""
function singlescaledbeta2D(betpar::ScaledBeta2DParams,xy::Array{<:Real,2},mcur::Vector{<:Real})

    # pre-allocate array of arrays (one for each y value)
    npts = size(xy,1)
    scbeta2d = zeros(Real,npts) #?? #zeros(eltype(mcur),npts)
    for l=1:npts
        ## xy[l,1] -> SDS concentration
        ## xy[l,2] -> protein concentration
        ycur = xy[l,2] # protein concentration
        xcur = xy[l,1]

        ## get the values of model parameters at y=ycur
        mode,kon,amp = getmodparbeta(betpar,mcur,ycur)

        # beta pdf along x
        ## xy[l,1] -> SDS concentration
        scbeta2d[l] = scaledbeta(mode,kon,betpar.a,betpar.b,amp,xcur)
    end

    return scbeta2d
end

###########################################################

"""
 $(TYPEDSIGNATURES)

Get the parameters (mode, confidence, amplitude) of the Beta functions for given y (protein concentration).
"""
function getmodparbeta(betpar::ScaledBeta2DParams,mcur::Vector{<:Real},ycur::Real)

    ##-------------------------------------------
    ## Linear: Straight line through 2 points:
    ##
    ## y = (y2-y1)/(x2-x1) * (x-x1) + y1
    ##

    ##-------------------------------------------
    ## Quadratic: to be implemented...
    ##
    ##
    ##

    ## limits of protein concentration
    x1 = betpar.ymin
    x2 = betpar.ymax
    @assert (x2-x1)!=0.0

    ## mode
    if betpar.modefuny=="linear"
        mode1,mode2 = mcur[betpar.idxmode:betpar.idxmode+1]
        mode = ((mode2-mode1)/(x2-x1))*(ycur-x1) + mode1
    elseif betpar.modefuny=="quadratic"
        error("Quadratic representation of mode along y not yet implemented.")
    end

    ## confidence parameter
    if betpar.konfuny=="constant"
        kon = mcur[betpar.idxkon]
    elseif betpar.konfuny=="linear"
        kon1,kon2 = mcur[betpar.idxkon:betpar.idxkon+1]
        kon = ((kon2-kon1)/(x2-x1))*(ycur-x1) + kon1
    elseif betpar.konfuny=="quadratic"
        error("Quadratic representation of confidence parameter along y not yet implemented.")
    end

    ## amplitude
    if betpar.ampfuny=="linear"
        amp1,amp2 = mcur[betpar.idxamp:betpar.idxamp+1]
        amp = ((amp2-amp1)/(x2-x1))*(ycur-x1) + amp1
    elseif betpar.ampfuny=="quadratic"
        error("Quadratic representation of amplitude parameter along y not yet implemented.")
    end

    return mode,kon,amp
end

###########################################################

"""
 $(TYPEDSIGNATURES)

1D modified scaled Beta function.
"""
@inline function scaledbeta(mo::Real,kon::Real,a::Real,b::Real,
                            amplscale::Real,x::Real)

    @assert b>a
    @assert a <= mo <= b
    ## keep kon at >2.0 instead of >=2.0 to avoid weird shapes...
    @assert kon>2.0

    rmo = (mo-a)/(b-a)
    α = rmo * (kon-2.0) + 1.0
    β = (1.0-rmo) * (kon-2.0) + 1.0
    @assert α>0
    @assert β>0

    betatyp = Beta(α,β)

    xresc = (x-a)/(b-a)
    tmpscbeval = pdf(betatyp,xresc)/(b-a)

    # btype = "absampl"

    # if btype == "pdfkind"
    #     ## pdf version with scaling, the value of amplscale does NOT correspond
    #     ##   to the actual value of the peak
    #     scbeval = amplscale * tmpscbeval

    # elseif btype == "absampl"
        ## here the value of amplscale corresponds to the actual value of the peak
        xrescmode = (mo-a)/(b-a)
        ## value at the mode
        peakval = pdf(betatyp,xrescmode)/(b-a)
        ## normalize first, then rescale
        scbeval = amplscale * (tmpscbeval/peakval)
    # end

    return scbeval
end

###############################################################################

###############################################################################

"""
$(TYPEDSIGNATURES)

Solve the inverse problem, i.e., fit the measured enthalpy data,
   using an Interior Point Newton method from the Optim.jl package.

# Arguments

- `betpar`: a struct containing the parameters for the Beta functions
            See [`ScaledBeta2DParams`](@ref)
- `dobs`: a struct containing the observed (measured) data and concentrations
            See [`ITCObsData`](@ref)
- `invCd`: inverse of the covariance matrix on observations (precision matrix)
- `mstart`: the starting model
- `lowconstr`: array of lower constraints for all parameters
- `upconstr`: array of upper constraints for all parameters
- `outdir`: output directory to save results
- `applynonlinconstr`=false: optional parameter determining whether to use or
                             not the nonlinear constraints
- `constrarea`=Inf: a positive real number defining the upper constraint for the value of area (enthalpy) at protein concentration equal to zero

# Returns

A structure holding Beta parameters and the solution in terms of mode, confidence parameter and amplitude.
It also saves all the setup of the problem and a set of parameters to an HDF5 file.

"""
function solveinvprob(betpar::ScaledBeta2DParams,dobs::ITCObsData,invCd::Matrix{<:Real},
                      mstart::Matrix{<:Real},lowconstr::Matrix{<:Real},upconstr::Matrix{<:Real},
                      outdir::String;
                      applynonlinconstr::Bool=false,constrarea::Real=Inf)

    # References

    # The algorithm was originally written by Tim Holy (@timholy, tim.holy@gmail.com).
    # J Nocedal, SJ Wright (2006), Numerical optimization, second edition. Springer.
    # A Wächter, LT Biegler (2006), On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Mathematical Programming 106 (1), 25-57.

    ### https://github.com/JuliaNLSolvers/NLSolversBase.jl/blob/0f07421dbef63e8f53ca658e6be2981c09ca9c5b/src/objective_types/constraints.jl
    ### Constraints
    #
    # Constraints are specified by the user as
    #    lx_i ≤   x[i]  ≤ ux_i  # variable (box) constraints
    #    lc_i ≤ c(x)[i] ≤ uc_i  # linear/nonlinear constraints
    # and become equality constraints with l_i = u_i. ±∞ are allowed for l
    # and u, in which case the relevant side(s) are unbounded.
    #
    # The user supplies functions to calculate c(x) and its derivatives.
    #
    # Of course we could unify the box-constraints into the
    # linear/nonlinear constraints, but that would force the user to
    # provide the variable-derivatives manually, which would be silly.

    # struct TwiceDifferentiableConstraints{F,J,H,T} <: Constraints
    #     c!::F # c!(storage, x) stores the value of the constraint-functions at x
    #     jacobian!::J # jacobian!(storage, x) stores the Jacobian of the constraint-functions
    #     h!::H   # h!(storage, x) stores the hessian of the constraint functions
    #     bounds::ConstraintBounds{T}
    # end

    @assert  size(mstart,1)==betpar.nummodpar

    ## Newton from Optim.jl
    #------------------------------------------------

    function closmisfitOPTIM(clmcur::Vector{<:Real})
        clmcur2d=reshape(clmcur,betpar.nummodpar,ncomp)
        msf = misfitfunctional(betpar,dobs,invCd,clmcur2d) #invCm,mprior,mcur)
        return msf
    end

    #------------------------------------------------

    function fun_grad!(gr,mcur)
        ## gradient of the misfit function
        gr .= ForwardDiff.gradient(closmisfitOPTIM,mcur)
        return nothing
    end

    #------------------------------------------------

    function fun_hess!(hess,mcur)
        ## Hessian of the misfit function
        hess .= ForwardDiff.hessian(closmisfitOPTIM,mcur)
        return nothing
    end

    #------------------------------------------------
    #------------------------------------------------

    if applynonlinconstr
        maxSDSc = 1e9

        function modeconstraint(betpar,ic::Integer,ccmcur::Vector{<:Real})
            ## compute the value of mode at prot. conc.=0
            ##   for a single beta function (ic)
            ycur = 0.0
            npar = betpar.nummodpar
            ncomp = div(length(ccmcur),npar)
            ccmcur2d = reshape(ccmcur,npar,ncomp)
            mode,kon,amp = getmodparbeta(betpar,ccmcur2d[:,ic],ycur)
            return mode
        end

        function konconstraint(betpar,ic::Integer,ccmcur::Vector{<:Real})
            ## compute the value of kon at prot. conc.=0
            ##   for a single beta function (ic)
            ycur = 0.0
            npar = betpar.nummodpar
            ncomp = div(length(ccmcur),npar)
            ccmcur2d = reshape(ccmcur,npar,ncomp)
            mode,kon,amp = getmodparbeta(betpar,ccmcur2d[:,ic],ycur)
            return kon
        end

        function ampconstraint(betpar,ic::Integer,ccmcur::Vector{<:Real})
            ## compute the value of amplitude at prot. conc.=0
            ##   for a single beta function (ic)
            ycur = 0.0
            npar = betpar.nummodpar
            ncomp = div(length(ccmcur),npar)
            ccmcur2d = reshape(ccmcur,npar,ncomp)
            mode,kon,amp = getmodparbeta(betpar,ccmcur2d[:,ic],ycur)
            return amp
        end

        function areaconstraint(betpar,ccmcur)
            ## compute the value of the area at prot. conc.=0
            ncomp = div(length(ccmcur),betpar.nummodpar)
            npar = betpar.nummodpar
            mcur2d = reshape(ccmcur,npar,ncomp)
            mar = misfareaenth(betpar,mcur2d)
            return mar
        end

        function angcoekonconstraint(betpar,ic::Integer,ccmcur)
            @assert betpar.konfuny=="linear"
            ## angular coefficient for kon parameter being negative
            ##    (shrinking functions), only if konfuny=="linear"
            npar = betpar.nummodpar
            ncomp = div(length(ccmcur),npar)
            ccmcur2d = reshape(ccmcur,npar,ncomp)
            konval = ccmcur2d[3:4,ic]
            ymin = betpar.ymin
            ymax = betpar.ymax
            kangcoe = (konval[2]-konval[1])/(ymax-ymin)
            return kangcoe
        end

        function modeatprotcon(betpar,ic::Integer,pconc::Real,ccmcur::Vector{<:Real})
            ## compute the value of mode at prot. conc.=0
            ##   for a single beta function (ic)
            ycur = pconc #0.0
            npar = betpar.nummodpar
            ncomp = div(length(ccmcur),npar)
            ccmcur2d = reshape(ccmcur,npar,ncomp)
            mode,kon,amp = getmodparbeta(betpar,ccmcur2d[:,ic],ycur)
            return mode
        end

        function nonintersectmodeconstr_below(betpar,ic::Integer,pconc::Real,ccmcur::Vector{<:Real})
            npar = betpar.nummodpar
            ncomp = div(length(ccmcur),npar)
            if ic==1
                # first Beta function
                modbelow = 0.0 # SDS=0.0 is the minimum that makes sense
            else
                # other Beta functions
                modbelow = modeatprotcon(betpar,ic-1,pconc,ccmcur) ## ic-1, Beta below
            end
            modcur = modeatprotcon(betpar,ic,pconc,ccmcur)
            ## difference in terms of SDS
            moddif = modcur-modbelow
            return moddif
        end

        function nonintersectmodeconstr_above(betpar,ic::Integer,pconc::Real,ccmcur::Vector{<:Real})
            npar = betpar.nummodpar
            ncomp = div(length(ccmcur),npar)
            if ic==ncomp
                # last Beta function
                modabove = maxSDSc #+Inf # SDS
            else
                # other Beta functions
                modabove = modeatprotcon(betpar,ic+1,pconc,ccmcur) ## ic-1, Beta below
            end
            modcur = modeatprotcon(betpar,ic,pconc,ccmcur)
            ## difference in terms of SDS
            moddif = modabove-modcur
            return moddif
        end

        function ampatpconc(betpar,ic::Integer,pconc,ccmcur::Vector{<:Real})
            ## compute the value of amplitude at prot. conc.=0
            ##   for a single beta function (ic)
            ycur = pconc
            npar = betpar.nummodpar
            ncomp = div(length(ccmcur),npar)
            ccmcur2d = reshape(ccmcur,npar,ncomp)
            mode,kon,amp = getmodparbeta(betpar,ccmcur2d[:,ic],ycur)
            return amp
         end

        function decramplconstraint(betpar,ic::Integer,ccmcur::Vector{<:Real})
            ##    trying to push towards zero amplitude at prot. conc.=0
            pc1=0.0
            pc2=betpar.ymax
            ## difference of amp between prot. conc. at betpar.ymax and 0.0
            amplat0 = ampatpconc(betpar,ic,pc1,ccmcur)
            amplatend = ampatpconc(betpar,ic,pc2,ccmcur)
            decamp = abs(amplatend)-abs(amplat0)
            # if amplat0<0.0
            #     #change sign
            #     decamp = -decamp
            # end
            # @show ic,decamp,amplat0,amplatend
            return decamp
        end

        #------------------------------------------------

        function funcon_c(betpar,ccmcur) #ccmcur::Vector{<:Real})
            ncomp = div(length(ccmcur),betpar.nummodpar)
            npar = betpar.nummodpar
            cstfun = Array{Any,1}(undef,0) #ncomp+ncomp+1)
            nonlincstnames = Array{Any,1}(undef,0)

            ## constraintS on mode at [protein]=0
            ycur = 0.0
            ifun = 1
            for ic=1:ncomp
                acfun = ccmcur->modeconstraint(betpar,ic,ccmcur)
                push!(cstfun,acfun)
                push!(nonlincstnames,"NONlinear constraint on mode at [protein]=0, Beta component #$ic")
                ifun+=1
            end

            ## constraints on kon at [protein]=0
            for ic=1:ncomp
                acfun = ccmcur->konconstraint(betpar,ic,ccmcur)
                push!(cstfun,acfun)
                push!(nonlincstnames,"NONlinear constraint on confidence at [protein]=0, Beta component #$ic")
                ifun+=1
            end

            # ## constraints on amp at [protein]=0
            # for ic=1:ncomp
            #     acfun = ccmcur->ampconstraint(betpar,ic,ccmcur)
            #     push!(cstfun,acfun)
            #     ifun+=1
            # end

            ## constraint on area at prot. conc.=0
            acfun = ccmcur-> areaconstraint(betpar,ccmcur)
            push!(cstfun,acfun)
            push!(nonlincstnames,"NONlinear constraint on value of area at prot. conc.=0")
            ifun+=1

            ## set constraints at prot. conc. 0.0 and ? such that
            ##  the mode lines cannot intersect
            for pconc in [0.0,betpar.ymax]
                ## constraints on the nonoverlapping mode lines BELOW
                for ic=1:ncomp
                    acfun = ccmcur->nonintersectmodeconstr_below(betpar,ic,pconc,ccmcur)
                    push!(cstfun,acfun)
                    push!(nonlincstnames,"NONlinear constraint on non-overlapping mode lines BELOW, Beta component #$ic")
                    ifun+=1
                end
                ## constraints on the nonoverlapping mode lines ABOVE
                for ic=1:ncomp
                    acfun = ccmcur->nonintersectmodeconstr_above(betpar,ic,pconc,ccmcur)
                    push!(cstfun,acfun)
                    push!(nonlincstnames,"NONlinear constraint on non-overlapping mode lines ABOVE, Beta component #$ic")
                    ifun+=1
                end
            end

            #############################################
            ## THE FOLLOWING MUST BE LAST CONSTRAINTS

            # angular coefficient for kon parameter being negative (shrinking functions), only if konfuny=="linear"
            if betpar.konfuny=="linear"
                for ic=1:ncomp
                    acfun = ccmcur->angcoekonconstraint(betpar,ic,ccmcur)
                    push!(cstfun,acfun)
                    push!(nonlincstnames,"NONlinear constraint on angular coefficient for confidence parameter, Beta component #$ic")
                    ifun+=1
                end
            end

            ## constraint on decreasing amplitude towards prot. con.=0
            if betpar.ampfuny=="linear"
               for ic=1:ncomp
                   acfun = ccmcur->decramplconstraint(betpar,ic,ccmcur)
                   push!(cstfun,acfun)
                   push!(nonlincstnames,"NONlinear constraint on decresing amplitude towards prot. conc.=0, Beta component #$ic")
                   ifun+=1
                end
            end

            return cstfun,nonlincstnames
        end

        ## actually create the functions
        cstfun,nonlincstnames = funcon_c(betpar,vec(mstart))

        #-------------------------------

        function con_c!(cst,ccmcur::Vector{<:Real})
            N=length(cstfun)
            for i=1:N
                cst[i] = cstfun[i](ccmcur)
                #@show i,cst[i]
            end
            return cst
        end

        ##===================================
        ## Set NONlinear constraints
        ncomp = size(mstart,2)
        #maxamp = 1.0

        lnlc = [zeros(ncomp)...,      # constr. on mode at [prot]==0
                2.5*ones(ncomp)...,   # constr. on kon at [prot]==0
                #-maxamp*ones(ncomp)..., # constr. on amp at [prot]==0
                0.0,           # constr. on area
                0.1.+zeros(ncomp)..., # constr. non-overlapping mode at prot. conc.=0 below
                0.1.+zeros(ncomp)..., # constr. non-overlapping mode at prot. conc.=0 above
                0.1.+zeros(ncomp)..., # constr. non-overlapping mode at prot. conc.=? below
                0.1.+zeros(ncomp)... ] # constr. non-overlapping mode at prot. conc.=? above

        unlc = [betpar.b*ones(ncomp)..., # constr. on mode at [prot]==0
                Inf*ones(ncomp)...,      # constr. on kon at [prot]==0
                #maxamp*ones(ncomp)...,  # constr. on amp at [prot]==0
                constrarea,              # constr. on area
                maxSDSc.*ones(ncomp)..., # constr. non-overlapping mode at prot. conc.=0 below
                maxSDSc.*ones(ncomp)..., # constr. non-overlapping mode at prot. conc.=0 above
                maxSDSc.*ones(ncomp)..., # constr. non-overlapping mode at prot. conc.=? below
                maxSDSc.*ones(ncomp)... ] # constr. non-overlapping mode at prot. conc.=? above

        if betpar.konfuny=="linear"
            # constr. on decreasing amplitude towards prot. con.=0
            lnlc = append!(lnlc,-Inf*ones(ncomp))
            # constr. on decreasing amplitude towards prot. con.=0
            unlc = append!(unlc,zeros(ncomp))
        end

        if betpar.ampfuny=="linear"
            # constr. on decreasing amplitude towards prot. con.=0
            lnlc = append!(lnlc,zeros(ncomp))
            # constr. on decreasing amplitude towards prot. con.=0
            unlc = append!(unlc,maxSDSc.*ones(ncomp))
        end

        ## init constraints
        cst = zeros(Real,length(lnlc))

        function closcon_c(mcur::Vector{<:Real})
            con_c!(cst,mcur)
            return cst
        end

        ##===================================

        function checknonlinconstr(betpar,cst,lnlc,unlc,mstvec,protein)
            cst = con_c!(cst,mstvec)
            fulfeachcon = lnlc.<=cst.<=unlc
            fulfillconstr = all(fulfeachcon)
            #     npar = betpar.nummodpar
            #     ncomp = div(length(mstvec),npar)
            #     mst = reshape(mstvec,npar,ncomp)
            #     betami = BetaMix2D(betpar,mst,protein)
            #     plotparamlines(betami)
            #     println(lnlc.<=cst.<=unlc)
            # end
            return fulfillconstr,fulfeachcon
        end

        ##===============================================
        #------------------------------------------------

        function con_jacob!(jac,mcur)
            #@show size(jac),size(mcur)
            jac .= ForwardDiff.jacobian(closcon_c,mcur)
            return nothing
        end

        #------------------------------------------------

        function con_hess!(hess,mcur, λ)
            ## IPNewton needs one Hessian per constraint
            # one Hessian for each single constraint
            N=length(λ)
            for i=1:N
                hess .+= λ[i] .* ForwardDiff.hessian(cstfun[i],mcur)
            end
            return nothing
        end

        #------------------------------------------------

    end # if applynonlinconstraints

    ##===================================================
    ## https://julianlsolvers.github.io/Optim.jl/stable/#examples/generated/ipnewton_basics/

    ## 2d -> 1D
    mstartvec = vec(mstart)
    ncomp = size(mstart,2)
    vlowconstr = vec(lowconstr)
    vupconstr = vec(upconstr)

    ##===================================================
    ## Some checks
    @assert issymmetric(invCd)
    @assert isposdef(invCd)
    @assert constrarea>=0.0

    mode,kon,amp = getmodparbeta(betpar,mstartvec,0.0)
    if mode<betpar.a
        betami = BetaMix2D(betpar,mstart,dobs.protein)
        plotparamlines(betami)
        error("solveinvprob(): mode<betpar.a at protein concentration=0.0 from the starting model.")
    end

    # check linear constraints
    fulfilleachLINcon = (lowconstr .<= mstart .<= upconstr)
    fulfillLINcon = all(fulfilleachLINcon)
    if !fulfillLINcon
        @warn("solveinvprob(): linear constraints not fulfilled by starting model")
        #println("Fulfilled linear constraints: ",string(fulfilleachLINcon))
        np,nc=size(fulfilleachLINcon)
        for ic=1:nc
            for ip=1:np
                if fulfilleachLINcon[ip,ic]==false
                    println("NOT FULFILLED: Linear constraint for parameter #$ip from Beta component #$ic")
                end
            end
        end
    end

    ##===================================================
    ## Objective function without constraints
    # df = TwiceDifferentiable(closmisfitOPTIM,mstartvec,autodiff=:forward)
    df = TwiceDifferentiable(closmisfitOPTIM,fun_grad!, fun_hess!, mstartvec)

    ## Add lin/nonlinear constraints
    if applynonlinconstr
        # check if starting model fulfills nonlinear constraints
        fullfillNONlinconstr,fulfeachNONlincon = checknonlinconstr(betpar,cst,lnlc,unlc,mstartvec,dobs.protein)
        if !fullfillNONlinconstr
            @warn("solveinvprob(): NONlinear constraints not fulfilled by starting model")
            #println("Fulfilled NONlinear constraints: ",string(fulfeachNONlincon))
            for i=1:length(fulfeachNONlincon)
                if fulfeachNONlincon[i]==false
                    println("##==> Cst #$i NOT FULFILLED: ",nonlincstnames[i])
                    println("##==> Cst #$i: lower constr. $(lnlc[i]) <= value: $(cst[i]) <= upper constr. $(unlc[i])")
                end
            end
        end
        # linear and NONlinear constraints
        dfc = TwiceDifferentiableConstraints(con_c!, con_jacob!, con_hess!,
                                             vlowconstr, vupconstr, lnlc, unlc)
    else
        # linear constraints only
        dfc = TwiceDifferentiableConstraints(vlowconstr, vupconstr)
    end

    ## Run the Newton inversion with box minimization
    println("\nRunning optimization with IPNewton...\n")
    # FAILURE using:  Optim.Options(store_trace=true,extended_trace=true) ????
    #result = optimize(df, dfc, mstartvec, IPNewton() ) #, Optim.Options(store_trace=true,extended_trace=true))
    # result = optimize(df, dfc, mstartvec, IPNewton() ) #, Optim.Options() )
    # to get Hessian: store_trace=true,extended_trace=true
    result = optimize(df, dfc, mstartvec, IPNewton(), Optim.Options(store_trace=true,extended_trace=true))

    rcon=Optim.converged(result)
    if rcon==false
        println("\n############################################")
        println(  "### FAILURE, convergence not attained!   ###")
        println(  "### Try with different input parameters. ###")
        println(  "############################################\n")
    end

    ##=========================================
    mpostvec = Optim.minimizer(result)
    ## https://julianlsolvers.github.io/Optim.jl/stable/#user/minimization/

    ## reshape mpost
    mpost = reshape(mpostvec,size(mstart))

    ## structure holding Beta parameters and the solution in terms of
    ##   mode, confidence parameter and amplitude
    betmix = BetaMix2D(betpar,mpost,dobs.protein)

    ## show some info
    println(result)

    ## inverse Hessian
    #@show size(result.trace[end].metadata["h(x)"])
    hessian = result.trace[end].metadata["h(x)"]

    ##================================
    ## save data
    try
        mkdir(outdir)
    catch
        nothing
    end

    ##-------------------------------------------------------
    # save in JLD2 format
    outfile1 = joinpath(outdir,dobs.protein*"_ITCinvresults.jld2")
    println("Saving results in JLD2 to $outfile1")
    jldopen(outfile1, "w") do fl
        fl["betamix"] = betmix
        fl["dobs"] = dobs
        fl["mstart"] = mstart
        fl["lowconstr"] = lowconstr
        fl["upconstr"] = upconstr
        fl["invCd"] = invCd
        fl["hessian"] = hessian
        if applynonlinconstr
            fl["lowNONlinconstr"] = lnlc
            fl["upNONlinconstr"] = unlc
            fl["NONlinconstrNAMES"] = nonlincstnames
        end
    end

    ##-------------------------------------------------------
    # save in HDF5 format
    # outfile1 = joinpath(outdir,dobs.protein*"_ITCresults.h5")
    # println("Saving results in HDF5 to $outfile1\n")

    # hf = h5open(outfile1,"w")
    # hf["protein"] = dobs.protein
    # hf["mpost"] = mpost
    # hf["enthalpy"] = dobs.enthalpy
    # hf["sdscon"] = dobs.sdsprotcon[:,1] # sdscon
    # hf["procon"] = dobs.sdsprotcon[:,2] # procon
    # hf["beta_a"] = betpar.a
    # hf["beta_b"] = betpar.b
    # hf["beta_ymin"] = betpar.ymax
    # hf["beta_ymax"] = betpar.ymin
    # hf["nummodpar"] = betpar.nummodpar
    # hf["mstart"] = mstart
    # hf["lowconstr"] = lowconstr
    # hf["upconstr"] = upconstr
    # hf["invCd"] = invCd
    # hf["modefuncy"] = betpar.modefuny
    # hf["konfuncy"]  = betpar.konfuny
    # hf["ampfuncy"] = betpar.ampfuny
    # close(hf)

    ##-------------------------------------------------------
    # save in plain text format
    outfile2 = joinpath(outdir,dobs.protein*"_ITCinvresults.dat")
    println("Saving results in text file to $outfile2\n")
    header = "# Output model parameters where each column represents a single Beta component for protein $(dobs.protein)"

    ofl = open(outfile2,"w")
    writedlm(ofl,[header])
    writedlm(ofl,mpost)
    close(ofl)

    return betmix
end

########################################################################
########################################################################

# function OLDsolveinvprob(betpar::ScaledBeta2DParams,dobs::ITCObsData,
#                       invCd::Matrix{<:Real},mstart::Matrix{<:Real},
#                       lowconstr::Vector{<:Real},upconstr::Vector{<:Real},
#                       outdir::String)

#     # References

#     # The algorithm was originally written by Tim Holy (@timholy, tim.holy@gmail.com).
#     # J Nocedal, SJ Wright (2006), Numerical optimization, second edition. Springer.
#     # A Wächter, LT Biegler (2006), On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Mathematical Programming 106 (1), 25-57.

#     ### https://github.com/JuliaNLSolvers/NLSolversBase.jl/blob/0f07421dbef63e8f53ca658e6be2981c09ca9c5b/src/objective_types/constraints.jl
#     ### Constraints
#     #
#     # Constraints are specified by the user as
#     #    lx_i ≤   x[i]  ≤ ux_i  # variable (box) constraints
#     #    lc_i ≤ c(x)[i] ≤ uc_i  # linear/nonlinear constraints
#     # and become equality constraints with l_i = u_i. ±∞ are allowed for l
#     # and u, in which case the relevant side(s) are unbounded.
#     #
#     # The user supplies functions to calculate c(x) and its derivatives.
#     #
#     # Of course we could unify the box-constraints into the
#     # linear/nonlinear constraints, but that would force the user to
#     # provide the variable-derivatives manually, which would be silly.

#     # struct TwiceDifferentiableConstraints{F,J,H,T} <: Constraints
#     #     c!::F # c!(storage, x) stores the value of the constraint-functions at x
#     #     jacobian!::J # jacobian!(storage, x) stores the Jacobian of the constraint-functions
#     #     h!::H   # h!(storage, x) stores the hessian of the constraint functions
#     #     bounds::ConstraintBounds{T}
#     # end

#     ## Newton from Optim.jl
#     #------------------------------------------------

#     function closmisfitOPTIM(clmcur::Vector{<:Real})
#         clmcur2d=reshape(clmcur,betpar.nummodpar,ncomp)
#         msf = misfitbeta2D(betpar,dobs,invCd,clmcur2d) #invCm,mprior,mcur)
#         return msf
#     end

#     #------------------------------------------------

#     function fun_grad!(gr,mcur)
#         ## gradient of the misfit function
#         gr .= ForwardDiff.gradient(closmisfitOPTIM,mcur)
#         return nothing
#     end

#     #------------------------------------------------

#     function fun_hess!(hess,mcur)
#         ## Hessian of the misfit function
#         hess .= ForwardDiff.hessian(closmisfitOPTIM,mcur)
#         return nothing
#     end

#     ##===================================================
#     ## Some checks
#     @assert issymmetric(invCd)
#     @assert isposdef(invCd)

#     ##===================================================
#     ## https://julianlsolvers.github.io/Optim.jl/stable/#examples/generated/ipnewton_basics/

#     ## 2d -> 1D
#     mstartvec = vec(mstart)
#     ncomp = size(mstart,2)

#     ## Objective function without constraints
#     df = TwiceDifferentiable(closmisfitOPTIM, fun_grad!, fun_hess!, mstartvec)

#     ## Add constraints
#     dfc = TwiceDifferentiableConstraints(lowconstr, upconstr)

#     ## Run the Newton inversion with box minimization
#     println("\nRunning optimization with IPNewton...")
#     result = optimize(df, dfc, mstartvec, IPNewton(), Optim.Options(store_trace=true))
#     mpostvec = Optim.minimizer(result)

#     ## reshape mpost
#     mpost = reshape(mpostvec,size(mstart))

#     ## structure holding Beta parameters and the solution in terms of
#     ##   mode, confidence parameter and amplitude
#     betmix = BetaMix2D(betpar,mpost,dobs.protein)

#     ## show some info
#     println(result)

#     ##================================
#     ## save data
#     try
#         mkdir(outdir)
#     catch
#         nothing
#     end

#     outfile = joinpath(outdir,dobs.protein*"_ITCresults.h5")
#     println("Saving results to $outfile\n")

#     hf = h5open(outfile,"w")
#     hf["protein"] = dobs.protein
#     hf["mpost"] = mpost
#     hf["enthalpy"] = dobs.enthalpy
#     hf["sdscon"] = dobs.sdsprotcon[:,1] # sdscon
#     hf["procon"] = dobs.sdsprotcon[:,2] # procon
#     hf["beta_a"] = betpar.a
#     hf["beta_b"] = betpar.b
#     hf["beta_ymin"] = betpar.ymax
#     hf["beta_ymax"] = betpar.ymin
#     hf["nummodpar"] = betpar.nummodpar
#     hf["mstart"] = mstart
#     hf["lowconstr"] = lowconstr
#     hf["upconstr"] = upconstr
#     hf["invCd"] = invCd
#     hf["modefuncy"] = betpar.modefuny
#     hf["konfuncy"]  = betpar.konfuny
#     hf["ampfuncy"] = betpar.ampfuny
#     close(hf)

#     return betmix
# end

# ####################################################################


text:
### Placing miscellaneuous utm functions here (helpers for common use cases)

# find the utm zone for a specified LLA position
# TODO: Make sure this is the same for all Ellipsoids
"""
Find the utm zone and band for the input type.  get_lat() and get_lon() must be defined for the input
"""
function utm_zone(lla)

    # int versions
    ilat = floor(Int64, bound_thetad(get_lat(lla)))
    ilon = floor(Int64, bound_thetad(get_lon(lla)))

    # get the latitude band
    band = max(-10, min(9,  fld((ilat + 80), 8) - 10))

    # and check for weird ones
    zone = fld((ilon + 186), 6);
    if ((band == 7) && (zone == 31) && (ilon >= 3))
        zone = 32
    elseif ((band == 9) && (ilon >= 0) && (ilon < 42))
        zone = 2 * fld((ilon + 183), 12) + 1
    end

    return (zone, band)

end

# TODO: overload with other ellipses as needed
"""
Find the SRID for utm zone that the lat lon point input lies in.  This is only valid for the WGS84 datum.

Returns the authority and the code
"""
utm_srid(lla::Union{LLA, LL}) = :EPSG, -1  # can't do it
function utm_srid(lla::Union{LLA{WGS84}, LL{WGS84}})

    # get the band and zone
    (zone, band) = utm_zone(lla)
    auth = :EPSG
    code = 32600 + (band > 0 ? 0 : 100) + zone
    return (auth, code)
end

ext: jl


text:
export bondstyle, nodestyle, textstyle

function bondstyle(shape::Symbol, properties...; r=0.075, dashed=false, dash_params=[1mm, 1mm])
    if !any(x->x isa Compose.Property{Compose.StrokePrimitive}, properties)
        properties = (properties..., stroke("black"))
    end
    if dashed
        properties = (properties..., strokedash(dash_params))
    end
    if shape == :arrow
        return compose(context(), arrow(), (context(), line()), properties...)
    end
    geometry = if shape == :line || shape==:default
        line()
    elseif shape == :rcurve
        curve((0.0, 0.0), (r, 0.0), (r, 1.0), (0.0,1.0))
    elseif shape == :lcurve
        curve((0.0, 0.0), (-r, 0.0), (-r, 1.0), (0.0,1.0))
    elseif shape == :dcurve
        curve((0.0, 0.0), (0.0, r), (1.0, r), (1.0,0.0))
    elseif shape == :ucurve
        curve((0.0, 0.0), (0.0, -r), (1.0, -r), (1.0,0.0))
    else
        error("shape $shape not defined.")
    end
    compose(context(), geometry, properties...)
end

function nodestyle(shape::Symbol, properties...; r=0.02, θ=0.0, w=1.0)
    geometry = if shape == :circle || shape==:default
        circle(0.0, 0.0, r)
    elseif shape == :diamond
        rot_ngon(θ, 0.0, 0.0, r, 4)
    elseif shape == :square
        rectangle(-r, -r, 2r, 2r)
    elseif shape == :rectangle
        rectangle(-r*w, -r, 2r*w, 2r)
    elseif shape == :triangle
        rot_ngon(θ, 0.0, 0.0, r, 3)
    else
        error("shape $shape not defined.")
    end
    compose(context(), geometry, properties...)
end

function textstyle(sym::Symbol, properties...; width=0.2, height=0.1)
    if sym == :default
        compose(context(), text(0.0, 0.0, "", hcenter, vcenter), properties...)
    elseif sym == :math
        compose(context(), mathjax(0.0, 0.0, width, height, ""), properties...)
    else
        error("style $sym not defined.")
    end
end

"""
    rot(a, b, θ)

rotate variables `a` and `b` by an angle `θ`
"""
function rot(a, b, θ)
    s, c = sincos(θ)
    a*c-b*s, a*s+b*c
end

function rot_polygon(θ::Real, points::AbstractVector)
    polygon(map(x->rot(x..., θ), points))
end

function rot_ngon(θ::Real, args...; kwargs...)
    c = first(ngon(args..., kwargs...).primitives)
    rot_polygon(θ, map(x->_value.(x), c.points))
end

_value(x::Compose.Add) = _value(x.a) + _value(x.b)
_value(x::Compose.Length) = x.value


text:
function osqp(P::AbstractArray{T},q::AbstractArray{U},
              A::AbstractArray{V},l::AbstractArray{W},u::AbstractArray{S},prob) where {T,U,V,W,S}

#      # TODO store and update this
#      prob = OSQP.Model()

     if  (T<:ForwardDiff.Dual || U<:ForwardDiff.Dual ||  V<:ForwardDiff.Dual || W<:ForwardDiff.Dual || S<:ForwardDiff.Dual)
         if (T<:ForwardDiff.Dual)
             Dtype = eltype(P)
             Ptype = typeof(P[1,1].partials)
         elseif (U<:ForwardDiff.Dual)
             Dtype = eltype(q)
             Ptype = typeof(q[1].partials)
         elseif (V<:ForwardDiff.Dual)
             Dtype = eltype(A)
             Ptype = typeof(A[1,1].partials)
         elseif (W<:ForwardDiff.Dual)
             Dtype = eltype(l)
             Ptype = typeof(l[1].partials)
         else
             Dtype = eltype(u)
             Ptype = typeof(u[1].partials)
         end

         if (T<:ForwardDiff.Dual)
             Pv = map(ForwardDiff.value, P)
             Pp = map(ForwardDiff.partials, P)
         else
             Pv = P
             Pp = zeros(Ptype, size(P))
         end
         if (U<:ForwardDiff.Dual)
             qv = map(ForwardDiff.value, q)
             qp = map(ForwardDiff.partials, q)
         else
             qv = q
             qp = zeros(Ptype, size(q))
         end
         if (V<:ForwardDiff.Dual)
             Av = map(ForwardDiff.value, A)
             Ap = map(ForwardDiff.partials, A)
         else
             Av = A
             Ap = zeros(Ptype, size(A))
         end
         if (W<:ForwardDiff.Dual)
             lv = map(ForwardDiff.value, l)
             lp = map(ForwardDiff.partials, l)
         else
             lv = l
             lp = zeros(Ptype, size(l))
         end
         if (S<:ForwardDiff.Dual)
             uv = map(ForwardDiff.value, u)
             up = map(ForwardDiff.partials, u)
         else
             uv = u
             up = zeros(Ptype, size(u))
         end

         m,n = size(Av)

         OSQP.setup!(prob; P=sparse(Pv), q=qv, A=sparse(Av), l=lv, u=uv, verbose=0)
#         OSQP.update!(prob, Px=sparse(Pv), q=qv, Ax=sparse(Av), l=lv, u=uv)
         results = OSQP.solve!(prob)

         xv = results.x
         xv[isnan.(xv)] .= 0.
         xv[isinf.(xv)] .= 0.
         λ_plus = max.(0., results.y)
         λ_minus = min.(0., results.y)
         λ_plus[isnan.(λ_plus)] .= 0.
         λ_plus[isinf.(λ_plus)] .= 0.
         λ_minus[isnan.(λ_minus)] .= 0.
         λ_minus[isinf.(λ_minus)] .= 0.

         Gv = vcat(Av, -Av)
         hv = vcat(uv, -lv)
         λv = vcat(λ_plus, -λ_minus)

         M = vcat(hcat(Pv, Gv'),
         hcat(Diagonal(λv)*Gv, Diagonal(Gv*xv - hv)))

         Apt = Ap[reshape(1:m*n,(m,n))']
         Gp = vcat(Ap, -Ap)
         Gpt = hcat(Apt, -Apt)
         hp = vcat(up, -lp)

         dres = vcat(Pp*xv + qp + Gpt*λv,
                     Diagonal(λv)*(Gp*xv - hp))

         xp = (-pinv(Matrix(M)) * dres)[1:n]

         # TODO you have to use the right type here...
         sol = map(Dtype, xv, xp)
     else
         OSQP.setup!(prob; P=sparse(P), q=q, A=sparse(A), l=l, u=u, verbose=0)
#          OSQP.update!(prob, Px=sparse(Pv), q=qv, Ax=sparse(Av), l=lv, u=uv)
         results = OSQP.solve!(prob)

         sol = results.x
     end

     return sol
end


text:
using Documenter, StanDiagnose

makedocs(
    modules = [StanDiagnose],
    format = Documenter.HTML(),
    checkdocs = :exports,
    sitename = "StanDiagnose.jl",
    pages = Any["index.md"]
)

deploydocs(
    repo = "github.com/goedman/StanDiagnose.jl.git",
)


text:
function Base.show(io::IO, p::Polynomial{T}) where {T}
    vars = variables(p)
    first = true
    exps = exponents(p)
    cfs = coefficients(p)

    m, n = size(exps)

    for i=1:n
        exp = exps[:, i]
        coeff = cfs[i]

        if (!first && show_plus(coeff))
            print(io, "+")
        end
        first = false

        if (coeff != 1 && coeff != -1) || exp == zeros(Int, m)
            show_coeff(io, coeff)
        elseif coeff == -1
            print(io, "-")
        end

        for (var, power) in zip(vars, exp)
            if power == 1
                print(io, pretty_var(var))
            elseif power > 1
                print(io, pretty_var(var), pretty_power(power))
            end
        end
    end

    if first
        print(io, zero(T))
    end
end

const SUBSCRIPT_TABLE = ('\u2080', '\u2081', '\u2082', '\u2083', '\u2084', '\u2085', '\u2086', '\u2087', '\u2088', '\u2089')
const SUPERSCRIPT_TABLE = ('\u2070', '\u00b9', '\u00b2', '\u00b3', '\u2074', '\u2075', '\u2076', '\u2077', '\u2078', '\u2079')

unicode_subscript(i::Int) = (SUBSCRIPT_TABLE[i + 1])
unicode_superscript(i::Int) = (SUPERSCRIPT_TABLE[i + 1])

function pretty_power(pow::Int)
    io_out = IOBuffer()
    _digits = digits(pow)
    for i in _digits
        print(io_out, unicode_superscript(i))
    end
    return reverse(String(take!(io_out)))
end

function pretty_var(var::String)
    m1 = match(r"([a-zA-Z]+)(?:_*)(\d+)", var)
    m2 = match(r"([a-zA-Z]+)(?:\[)(\d+)(?:\])", var)
    if isnothing(m1) && isnothing(m2)
        return var
    end
    m = isnothing(m1) ? m2 : m1
    io_out = IOBuffer()
    print(io_out, m.captures[1])
    index = parse(Int, m.captures[2])
    _digits = reverse(digits(index))
    for i in _digits
        print(io_out, unicode_subscript(i))
    end
    return String(take!(io_out))
end
pretty_var(var) = pretty_var(string(var))

# helpers
show_plus(x::Real) = x >= 0
show_plus(x::Complex) = x != -1

show_coeff(io::IO, x::Real) = print(io, x)
function show_coeff(io::IO, x::Complex)
    if imag(x) == 0.0
        print(io, convert(Float64, x))
    else
        print(io, "($(x))")
    end
end


text:
using Documenter, MatLang

makedocs(;
    modules=[MatLang],
	format = Documenter.HTML(
        prettyurls = prettyurls = get(ENV, "CI", nothing) == "true",
    ),
    pages=[
        "Home" => "index.md",
		"Functions" =>
			["Functions List" => "functions.md",
			"Functions Explanations and Examples" =>
				["Entering Commands" => "functions/Entering_Commands.md",
				"Matrices and Arrays" => "functions/Matrices_and_Arrays.md",
				"Numeric Types" => "functions/Numeric_Types.md"],],
        "Development - Contribution" => "development.md",
        "Native Julia noteworthy differences from MATLAB" => "juliavsmatlab.md"
    ],
    repo="https://github.com/juliamatlab/MatLang/blob/{commit}{path}#L{line}",
    sitename="MatLang",
    authors="Amin Yahyaabadi - JuliaMatlab",
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(;
    repo="github.com/juliamatlab/MatLang.git",
)


text:
module Input

  import MPI
  using ODLCommonTools
  using PETSc2
  using Utils
  export read_input, make_input, read_input_file, registerOptionsChecker,
         printOpts

  include("physics_specific.jl")
  include("read_input.jl")
  include("make_input.jl")
  include("known_keys.jl")
  include("utils.jl")

end  # end module


text:
struct ImmutablePipeline{N,T<:Tuple} <: Pipeline
    operations::T

    function ImmutablePipeline{N}(ops::NTuple{N,Operation}) where N
        new{N,typeof(ops)}(ops)
    end
end
ImmutablePipeline() = throw(MethodError(ImmutablePipeline, ()))
ImmutablePipeline(::Tuple{}) = throw(MethodError(ImmutablePipeline, ((),)))
ImmutablePipeline(ops::NTuple{N,Operation}) where {N} = ImmutablePipeline{N}(ops)
ImmutablePipeline(ops::Vararg{Operation,N}) where {N} = ImmutablePipeline{N}(ops)
(::Type{Pipeline})(ops::NTuple{N,Operation}) where {N} = ImmutablePipeline(ops)
(::Type{Pipeline})(ops::Vararg{Operation,N}) where {N} = ImmutablePipeline(ops)

@inline Base.length(p::ImmutablePipeline{N}) where {N} = N
@inline operations(p::ImmutablePipeline) = p.operations
@inline operations(tup::NTuple{N,Operation}) where {N} = tup

# Allow specifying pipelines using the "op1 |> op2 |> op3" syntax.
Base.:(|>)(op1::Operation, op2::Operation) =
    ImmutablePipeline(op1, op2)
Base.:(|>)(p1::ImmutablePipeline, op2::Operation) =
    ImmutablePipeline(operations(p1)..., op2)
Base.:(|>)(op1::Operation, p2::ImmutablePipeline) =
    ImmutablePipeline(op1, operations(p2)...)
Base.:(|>)(p1::ImmutablePipeline, p2::ImmutablePipeline) =
    ImmutablePipeline(operations(p1)..., operations(p2)...)
Base.:(|>)(op1::Operation, buffer::AbstractArray) =
    ImmutablePipeline(op1, CacheIntermediate(buffer))
Base.:(|>)(p1::ImmutablePipeline, buffer::AbstractArray) =
    ImmutablePipeline(operations(p1)..., CacheIntermediate(buffer))

function Base.show(io::IO, pipeline::Pipeline)
    n = length(pipeline)
    ops = operations(pipeline)
    if get(io, :compact, false)
        for (i, op) in enumerate(ops)
            showconstruction(io, op)
            i < n && print(io, " |> ")
        end
    else
        k = length("$n")
        print(io, "$n-step $(typeof(pipeline).name):")
        for (i, op) in enumerate(ops)
            println(io)
            print(io, lpad(string(i), k+1, " "), ".) ")
            show(IOContext(io, :compact => true), op)
        end
    end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

# ToDo: Add literature references to AdaptiveMHTuning docstring.

"""
    struct AdaptiveMHTuning <: MHProposalDistTuning

Adaptive MCMC tuning strategy for Metropolis-Hastings samplers.

Adapts the proposal function based on the acceptance ratio and covariance
of the previous samples.

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)
"""
@with_kw struct AdaptiveMHTuning <: MHProposalDistTuning
    "Controls the weight given to new covariance information in adapting the
    proposal distribution."
    λ::Float64 = 0.5

    "Metropolis-Hastings acceptance ratio target, tuning will try to adapt
    the proposal distribution to bring the acceptance ratio inside this interval."
    α::IntervalSets.ClosedInterval{Float64} = ClosedInterval(0.15, 0.35)

    "Controls how much the spread of the proposal distribution is
    widened/narrowed depending on the current MH acceptance ratio."
    β::Float64 = 1.5

    "Interval for allowed scale/spread of the proposal distribution."
    c::IntervalSets.ClosedInterval{Float64} = ClosedInterval(1e-4, 1e2)

    "Reweighting factor. Take accumulated sample statistics of previous
    tuning cycles into account with a relative weight of `r`. Set to
    `0` to completely reset sample statistics between each tuning cycle."
    r::Real = 0.5
end

export AdaptiveMHTuning

mutable struct ProposalCovTuner{
    S<:MCMCBasicStats
} <: AbstractMCMCTunerInstance
    config::AdaptiveMHTuning
    stats::S
    iteration::Int
    scale::Float64
end

(tuning::AdaptiveMHTuning)(chain::MHIterator) = ProposalCovTuner(tuning, chain)

function ProposalCovTuner(tuning::AdaptiveMHTuning, chain::MHIterator)
    m = totalndof(getdensity(chain))
    scale = 2.38^2 / m
    ProposalCovTuner(tuning, MCMCBasicStats(chain), 1, scale)
end

_approx_cov(target::Distribution) = cov(unshaped(target))
_approx_cov(target::DistLikeDensity) = cov(target)
_approx_cov(target::AbstractPosteriorDensity) = cov(getprior(target))
_approx_cov(target::BAT.TransformedDensity{<:Any,<:BAT.DistributionTransform}) =
    BAT._approx_cov(target.trafo.target_dist)

function tuning_init!(tuner::ProposalCovTuner, chain::MHIterator)
    Σ_unscaled = _approx_cov(getdensity(chain))
    Σ = Σ_unscaled * tuner.scale

    next_cycle!(chain) # ToDo: This would be better placed in the burn-in algorithm
    chain.proposaldist = set_cov(chain.proposaldist, Σ)

    nothing
end

function tuning_postinit!(tuner::ProposalCovTuner, chain::MHIterator, samples::DensitySampleVector)
    # The very first samples of a chain can be very valuable to init tuner
    # stats, especially if the chain gets stuck early after:
    stats = tuner.stats
    append!(stats, samples)
end

function tuning_update!(tuner::ProposalCovTuner, chain::MHIterator, samples::DensitySampleVector)
    stats = tuner.stats
    stats_reweight_factor = tuner.config.r
    reweight_relative!(stats, stats_reweight_factor)
    # empty!.(stats)
    append!(stats, samples)

    config = tuner.config

    α_min = minimum(config.α)
    α_max = maximum(config.α)

    c_min = minimum(config.c)
    c_max = maximum(config.c)

    β = config.β

    t = tuner.iteration
    λ = config.λ
    c = tuner.scale
    Σ_old = Matrix(get_cov(chain.proposaldist))

    S = convert(Array, stats.param_stats.cov)
    a_t = 1 / t^λ
    new_Σ_unscal = (1 - a_t) * (Σ_old/c) + a_t * S

    α = eff_acceptance_ratio(chain)

    max_log_posterior = stats.logtf_stats.maximum

    if α_min <= α <= α_max
        chain.info = MCMCIteratorInfo(chain.info, tuned = true)
        @debug "MCMC chain $(chain.info.id) tuned, acceptance ratio = $(Float32(α)), proposal scale = $(Float32(c)), max. log posterior = $(Float32(max_log_posterior))"
    else
        chain.info = MCMCIteratorInfo(chain.info, tuned = false)
        @debug "MCMC chain $(chain.info.id) *not* tuned, acceptance ratio = $(Float32(α)), proposal scale = $(Float32(c)), max. log posterior = $(Float32(max_log_posterior))"

        if α > α_max && c < c_max
            tuner.scale = c * β
        elseif α < α_min && c > c_min
            tuner.scale = c / β
        end
    end

    Σ_new = new_Σ_unscal * tuner.scale

    chain.proposaldist = set_cov(chain.proposaldist, Σ_new)
    tuner.iteration += 1

    nothing
end


text:
function add!(n::NetworkS,extra::Extracellular)

    extra.v=zeros(Float64,length(n.t))
    coeffs=Array(Extra_coeffs,0)

    for j=1:length(fieldnames(n.neur)),k in getfield(n.neur,j)
        mycoeffs=extracellular(extra,k,0.3)
        push!(coeffs,mycoeffs)
    end

    add_extra_(n,extra,coeffs)
end

add_extra_(n::Network,e::Extracellular,c::Array{Extra_coeffs,1})=(e=typeof(e)(e.xyz,c,e.v);push!(n.extra,e))

add!(n::NetworkS,s::Stim)=add_stim_(n,s)

function add_stim_(n::Network,s::Stim)
    myis=zeros(Float64,length(n.t))

    startind=findfirst(n.t.>s.tstart)
    endind=findfirst(n.t.>s.tstop)

    myis[startind:endind]=s.Is[1]

    s.Is=myis

    push!(n.stim,s)
end

add!(n::NetworkS,intra::Intracellular)=add_intra_(n,intra)

add_intra_(n::Network,i::Intracellular)=(i.v=zeros(Float64,length(n.t));push!(n.intra,i))

function run!(n::NetworkS,init=false)

    #get initial conditions if uninitialized
    if init==true
        init!(n)
    end

    for i=1:length(n.t)

        for j=1:length(n.stim)
            @inbounds getfield(n.neur,n.stim[j].mtype)[n.stim[j].neur].rhs[n.stim[j].node]+=n.stim[j].Is[i]
        end

        if n.helper.flags[1]
            count=1
            for j=1 : length(fieldnames(n.neur))
                for k=1:length(getfield(n.neur,j))
                    @inbounds main(getfield(n.neur,j)[k])
                    for l=1:length(n.extra)
                        @inbounds n.extra[l].v[i]+=a_mult_b(n.extra[l].coeffs[count],getfield(n.neur,j)[k])
                    end
                    count+=1
                end
            end
        else
            @inbounds main(n.neur)
        end

        for j=1:length(n.intra)
            @inbounds n.intra[j].v[i]=getfield(n.neur,n.intra[j].mtype)[n.intra[j].neur].v[n.intra[j].node]
        end
    end

    nothing
end

function init!(n::NetworkS)

    @inbounds initialcon!(n.neur)

    #set up flags
    if length(n.extra)==0
        n.helper.flags[1]=false
    else
        n.helper.flags[1]=true
    end

    if length(n.intra)==0
        n.helper.flags[2]=false
    else
        n.helper.flags[2]=true
    end

    if length(n.stim)==0
        n.helper.flags[3]=false
    else
        n.helper.flags[3]=true
    end
end

function a_mult_b(x::Extra_coeffs,n::Neuron)
    c=0.0
    for i=1:length(x.c)
        @inbounds c+=x.c[i]*n.i_vm[x.inds[i]]
    end
    c
end

max_stars_repo_name: paulmthompson/JNeuron

max_issues_repo_name: paulmthompson/JNeuron

max_forks_repo_name: paulmthompson/JNeuron


text:
# Logistic Regression with Gradient Descent applied to MNIST Dataset

push!(LOAD_PATH, ".")
push!(LOAD_PATH, "../common")

using MNIST
importall logistic_regression
importall optimization

# parameters
alpha = 0.5
max_iter = 20
min_err = 0.0001

# load MNIST data
println("Loading MNIST dataset...")
X_train, y_train = traindata()
X_test, y_test = testdata()

# filter MNIST data, since it is binary classification
filt_train, filt_test =  y_train .< 2.0, y_test .< 2.0
X_train, y_train = X_train[:,filt_train], y_train[filt_train]
X_test, y_test = X_test[:,filt_test], y_test[filt_test]

# add ones as bias to feature vectors, feature scaling
X_train = [255*ones(1, size(X_train, 2)); X_train] / 255
X_test = [255*ones(1, size(X_test, 2)); X_test] / 255

# transpose results arrays
y_train, y_test = y_train', y_test'

# theta initialization
theta = init_theta(size(X_train, 1))

# training
println("Training...")
tic()
theta, history = gradient_descent(h, J, grad, X_train, y_train, theta,
                                  alpha, max_iter, min_err)
toc()

println("train accuracy: ", accuracy(theta, X_train, y_train))
println("test accuracy: ", accuracy(theta, X_test, y_test))

max_stars_repo_name: ilkerkesen/ufldl-tutorial

max_issues_repo_name: ilkerkesen/ufldl-tutorial

max_forks_repo_name: ilkerkesen/ufldl-tutorial


text:
using Test
using CloudCovErr

tests = [
    "tst_preprocess.jl",
    "tst_cov.jl",
    "tst_perstar.jl",
    "tst_decam.jl"
]

for t in tests
    @testset "$t" begin
        include(t)
    end
end


text:
#= ForwardDiff is removed from dependencies in this version.
import MechanicalSketch
import MechanicalSketch: color_with_lumin, background, O, WI, HE, EM, FS, finish,
       PALETTE, color_from_palette, setfont, settext, empty_figure, sethue, @layer
import MechanicalSketch: dimension_aligned
import MechanicalSketch: ComplexQuantity, generate_complex_potential_source, generate_complex_potential_vortex
import MechanicalSketch: @import_expand, BinLegend, BinLegendVector, draw_legend
import MechanicalSketch: place_image, set_scale_sketch, lenient_min_max
import MechanicalSketch: ∙, ∇_rectangle, x_y_iterators_at_pixels

let
BACKCOLOR = color_with_lumin(PALETTE[8], 10)
function restart()
    empty_figure(filename = joinpath(@__DIR__, "test_23.png"))
    background(BACKCOLOR)
    sethue(PALETTE[8])
end
restart()

if !@isdefined m²
    @import_expand(m, s)
end

"Source position"
p_source = complex(3.0, 0.0)m
"Flow rate, 2d flow"
q_source = 1.0m²/s
"Vortex position"
p_vortex = complex(0.0, 1.0)m
"Vorticity, 2d flow"
K = 1.0m²/s / 2π

ϕ_vortex =  generate_complex_potential_vortex(; pos = p_vortex, vorticity = K)
ϕ_source = generate_complex_potential_source(; pos = p_source, massflowout = q_source)
ϕ_sink = generate_complex_potential_source(; pos = -p_source, massflowout = -q_source)

"""
    ϕ(p::ComplexQuantity)
    → Quantity
2d velocity potential function. Complex quantity domain, real quantity range.
"""
ϕ(p) = ϕ_vortex(p) + ϕ_source(p) + ϕ_sink(p)

ϕ(p_source)
ϕ(p_vortex + p_source)

# Plot the real-valued function
OT = O + (-EM, -0.25HE + EM)
physwidth = 10.0m
physheight = 0.4 * physwidth
set_scale_sketch(physwidth, round(Int, WI * 2 / 3))
xs, ys = x_y_iterators_at_pixels(;physwidth, physheight)
A = [ϕ(complex(x, y)) for y in ys, x in xs]
toplegend = BinLegend(;maxlegend = maximum(A), minlegend = -maximum(A),
                      noofbins = 256, name = :Potential)
colormat = toplegend.(A);
upleftpoint, lowrightpoint = place_image(OT, colormat)
legendpos = lowrightpoint + (EM, 0) + (0.0m, physheight)
draw_legend(legendpos, toplegend)

# Text
setfont("DejaVu Sans", FS)
str = "ϕ: Z ↣ R  is the flow potenial"
settext(str, O + (-WI/2 + EM, -0.5HE + 2EM), markup = true)
settext("Source", OT + (real(p_source), imag(p_source)))
settext("Vortex", OT + (real(p_vortex), imag(p_vortex)))
@layer begin
    sethue(BACKCOLOR)
    settext("Sink", OT - (real(p_source), imag(p_source)))
end
str = "∇ϕ: Z ↣ Z  is the flow gradient, aka velocity vectors"
settext(str, O + (-WI/2 + EM, 1.5EM), markup = true)

# Plot the complex-valued function, aka velocity vectors.
B = begin
        unclamped = ∇_rectangle(ϕ; physwidth, physheight)
        map(unclamped) do u
            hypot(u) > 0.5m/s ? NaN∙u : u
        end
    end

OB = O + (-EM, + 0.25HE + 0.5EM )
mi, ma = lenient_min_max(B)
botlegend = BinLegendVector(;operand_example = first(B),
        max_magn_legend = ma, noof_magn_bins = 30, noof_ang_bins = 36,
        name = :Velocity)
colormat = botlegend.(B)
ulp, lrp = place_image(OB, colormat)
legendpos = lrp + (EM, 0) + (0.0m, physheight)
draw_legend(legendpos, botlegend)

# Position indicator
dimension_aligned(OB + (-physwidth / 2, physheight / 2), OB + (physwidth / 2, physheight / 2))
dimension_aligned(OB + p_vortex, OB)
dimension_aligned(OB + (-physwidth / 2, - physheight / 2 ),  OB +  (-physwidth / 2, physheight / 2 ))

finish()
end #let
=#


text:
# This file was generated, do not modify it. # hide
figure(figsize=(8, 6))
hist(data.Wind_gen, color="blue", edgecolor = "white", bins=50,
     density=true, alpha=0.5)
xlabel("Wind power generation (MWh)", fontsize=14)
ylabel("Frequency", fontsize=14)

savefig(joinpath(@OUTPUT, "hist_wind.svg")) # hide


text:
const Cube = Array{T, 3} where T

mutable struct State{F <: AbstractFloat}
  Z::Matrix{Bool} # Dim: J x K. Z[j,k] ∈ {0, 1}. true => 1, false => 0
  delta::Dict{Bool, Vector{F}} # delta_z dim: L[z]
  alpha::F
  v::Vector{F} # dim: K
  W::Matrix{F} # dim: I x K
  sig2::Vector{F} # dim: I
  eta::Dict{Bool, Cube{F}} # eta_zik dim: L[z]
  lam::Vector{Vector{Int8}} # Array of Array. lam[1:I] ∈ {1,...,K}
  gam::Vector{Matrix{Int8}} # gam_ij dim: N[i]
  y_imputed::Vector{Matrix{F}} # y_ij dim: N[i]
  eps::Vector{F} # dim I
  State{F}() where {F <: AbstractFloat} = new()
end

#= Note:
genInitialState(c::Constants, d::Data)
is in Constants.jl.
=#

function compress(state::State)
  warn_msg = "WARNING: The `compress` function has not been fully tested "
  warn_msg *= "and may result in errors!"
  println(warn_msg)

  if typeof(state) == State{Float32}
    return state
  else
    new_state = State{Float32}()
    new_state.Z = Matrix{Bool}(state.Z)
    new_state.delta = Dict{Bool, Vector{Float32}}(state.delta)
    new_state.alpha = Float32(state.alpha)
    new_state.v = Float32.(state.v)
    new_state.W = Float32.(state.W)
    new_state.sig2 = Float32.(state.sig2)
    new_state.eta = Dict{Bool, Cube{Float32}}(state.eta)
    new_state.lam = Vector{Vector{Int8}}(state.lam)
    new_state.gam = Vector{Matrix{Int8}}(state.gam)
    new_state.eps = Float32.(state.eps)
    new_state.y_imputed = Vector{Matrix{Float32}}(state.y_imputed)
    return new_state
  end
end


text:
export write_data

function write_data(filename::String, args...)
    error("Please check function signature. \n $(methods(write_data)))")
end

function write_data_peridigm(filename::String, x::Matrix, type::Vector,  vol::Vector)
    file = open(filename, "w+")
    N = size(x, 2)
    write(file, "$N \n\n")
    for j in 1:size(x, 2)
        t = type[j]
        v = vol[j]
        a,b,c = x[1,j],x[2,j],x[3,j]
        write(file, "$a, $b, $c, $t, $v ")
        write(file,"\n")
    end
    close(file)
end

function write_data(filename::String, x::Matrix, v::Matrix, type::Vector,  vol::Vector)
    file = open(filename, "w+")
    N = size(x, 2)
    write(file, "$N \n\n")
    for j in 1:size(x, 2)
        t = type[j]
        vol_ = vol[j]
        a,b,c = x[:, j]
        d,e,f = v[:, j]
        write(file, "$t, $a, $b, $c, $d, $e, $f, $vol_")
        write(file,"\n")
    end
    close(file)
end

max_stars_repo_name: ravinderbhattoo/PDMesh

max_issues_repo_name: ravinderbhattoo/PDMesh

max_forks_repo_name: ravinderbhattoo/PDMesh