##
# Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
#
# File :      portfolio_5_card.jl
#
# Description :  Implements a basic portfolio optimization model.

using Mosek
function portfolio( mu :: Vector{Float64},
                    x0 :: Vector{Float64},
                    w  :: Float64,
                    gamma :: Float64,
                    GT :: Array{Float64,2},
                    K  :: Int)
    (k,n) = size(GT)
    maketask() do task
        # Use remote server: putoptserverhost(task,"http://solve.mosek.com:30080")
        # Directs the log task stream
        # putstreamfunc(task,MSK_STREAM_LOG,msg -> print(msg))

        totalBudget = sum(x0)+w

        #Offset of variables into the API variable.
        x_ofs = 0
        y_ofs = n
        z_ofs = 2*n

        # Constraints offsets

        # Holding variable x of length n
        # No other auxiliary variables are needed in this formulation
        appendvars(task,3*n)

        # Setting up variable x
        for j in 1:n
            # Optionally we can give the variables names
            putvarname(task, x_ofs+j, "x[$(j)]");
            putvarname(task, y_ofs+j, "y[$(j)]");
            putvarname(task, z_ofs+j, "z[$(j)]");
            putvartype(task, y_ofs+j, MSK_VAR_TYPE_INT);
        end
        # No short-selling - x^l = 0, x^u = inf
        putvarboundsliceconst(task,x_ofs+1,x_ofs+n+1, MSK_BK_LO, 0.0, Inf);
        putvarboundsliceconst(task,y_ofs+1,y_ofs+n+1, MSK_BK_RA, 0.0, 1.0);
        putvarboundsliceconst(task,z_ofs+1,z_ofs+n+1, MSK_BK_FR, -Inf, Inf);

        # One linear constraint: total budget
        let coni = getnumcon(task)
            appendcons(task,1);
            putconname(task,coni+1,"budget")

            consub = Int32[ coni+1 for i in 1:n ]
            putaijlist(task,consub,[x_ofs+1:x_ofs+n...], ones(n))
            putconbound(task, coni+1, MSK_BK_FX, totalBudget, totalBudget)
        end

        let coni = getnumcon(task)
            appendcons(task,1)
            putarow(task,coni+1,[y_ofs+1:y_ofs+1+n...],ones(n));
            putconbound(task,coni+1,MSK_BK_UP, 0.0, K)
            putconname(task,coni+1,"cardinality")
        end

        let coni = getnumcon(task)
            appendcons(task,2*n)
            for i in 1:n
                putconname(task,coni+i,"zabs1[$i]")
                putconname(task,coni+n+i,"zabs2[$i]")

                putaij(task,coni + i, x_ofs + i, -1.0)
                putaij(task,coni + i, z_ofs + i, 1.0)
                putconbound(task,coni + i, MSK_BK_LO, -x0[i], Inf)

                putaij(task,coni + n + i, x_ofs + i, 1.0)
                putaij(task,coni + n + i, z_ofs + i, 1.0)
                putconbound(task,coni + n + i, MSK_BK_LO, x0[i], Inf)
            end
        end

        # - Switch
        let coni = getnumcon(task)
            appendcons(task,n)
            for i in 1:n
                putconname(task,coni+i, "switch[$i]")
                putaij(task,coni + i, z_ofs + i, 1.0)
                putaij(task,coni + i, y_ofs + i, -totalBudget)
                putconbound(task,coni + i, MSK_BK_UP, -Inf, 0.0)
            end
        end

        let afei = getnumafe(task),
            acci = getnumacc(task)
            # Input (gamma, GTx) in the AFE (affine expression) storage
            # We need k+1 rows
            appendafes(task,k + 1)
            # The first affine expression = gamma
            putafeg(task,afei+1, gamma)
            # The remaining k expressions comprise GT*x, we add them row by row
            # In more realisic scenarios it would be better to extract nonzeros and input in sparse form


            subj = [1:n...]
            for i in 1:k
                putafefrow(task,afei+i+1, subj, GT[i,:])
            end

            # Input the affine conic constraint (gamma, GT*x) \in QCone
            # Add the quadratic domain of dimension k+1
            qdom = appendquadraticconedomain(task,k + 1)
            # Add the constraint
            appendaccseq(task,qdom,1,nothing)
            putaccname(task,acci+1, "risk")
        end


        # Objective: maximize expected return mu^T x
        putclist(task,[x_ofs+1:x_ofs+n...],mu)
        putobjsense(task,MSK_OBJECTIVE_SENSE_MAXIMIZE)

        optimize(task)

        if getsolsta(task, MSK_SOL_ITG) != MSK_SOL_STA_INTEGER_OPTIMAL
            # See https://docs.mosek.com/latest/juliaapi/accessing-solution.html about handling solution statuses.
            error("Solution not optimal")
        end
        
        # Display solution summary for quick inspection of results
        solutionsummary(task,MSK_STREAM_LOG)

        writedata(task,"portfolio_5_card-$K.ptf");

        # Read the results
        xx = getxxslice(task,MSK_SOL_ITG, x_ofs+1,x_ofs+n+1)
        expret = mu' * xx

        (xx,expret)
    end
end # portfolio()

let w = 1.0,
    mu = [0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379],
    x0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
    GT = [ 0.30758 0.12146 0.11341 0.11327 0.17625 0.11973 0.10435 0.10638
           0.      0.25042 0.09946 0.09164 0.06692 0.08706 0.09173 0.08506
           0.      0.      0.19914 0.05867 0.06453 0.07367 0.06468 0.01914
           0.      0.      0.      0.20876 0.04933 0.03651 0.09381 0.07742
           0.      0.      0.      0.      0.36096 0.12574 0.10157 0.0571
           0.      0.      0.      0.      0.      0.21552 0.05663 0.06187
           0.      0.      0.      0.      0.      0.      0.22514 0.03327
           0.      0.      0.      0.      0.      0.      0.      0.2202 ],
    (k,n) = size(GT),
    gamma  = 0.25

    for K in 1:n
        res = Float64[]
        let (xx,expret) = portfolio(mu,x0,w,gamma,GT,K)
            println("cardinality = $K")
            println("\tExpected return $(expret) for gamma $(gamma)")
            println("\tSolution vector = $(xx)")
        end
    end
end
