using JuMP, PATHSolver#, NamedArrays, DataFrames
 mcp = Model(PATHSolver.Optimizer)

    @variable(mcp, X, (start = 1))
    @variables(mcp, begin
        Y, (start = 1)
        W, (start = 1)
        PX, (start = 1)
        PY, (start = 1)
        PW, (start = 1)
        PL, (start = 1)
        PK, (start = 1)
        CONS, (start = 200)
        tax_pl in JuMP.Parameter(0)
        tax_pk in JuMP.Parameter(0)
    end)

    cost_X_va = (PL*(1+tax_pl))^(40/100) * (PK*(1+tax_pk))^(60/100)
    cost_X = (20/120 * PY^(1-.5) + 100/120 * cost_X_va^(1-.5))^(1/(1-.5))
    revenue_X = PX

    cost_Y_va = (PL)^(60/100) * (PK)^(40/100)
    cost_Y = (20/120 * PX^(1-.75) + 100/120 * cost_Y_va^(1-.75))^(1/(1-.75))
    revenue_Y = PY

    cost_W = PY^.5*PX^.5
    revenue_W = PW

        # Zero Profit
        @constraints(mcp, begin
            zp_X, 120*cost_X - 120*revenue_X ⟂ X
            zp_Y, 120*cost_Y - 120*revenue_Y ⟂ Y
            zp_W, 200*cost_W - 200*revenue_W ⟂ W
        end)

        # Market Clearance
        @constraints(mcp, begin
            mc_PX, 120*X - 20*Y*(cost_Y/PX)^.75 - 100*W*(cost_W/PX) ⟂ PX
            mc_PY, -20*X*(cost_X/PY)^.5 + 120*Y - 100*W*(cost_W/PY) ⟂ PY
            mc_PW, 200*W - CONS/W ⟂ PW
            mc_PL, -40*X*(cost_X/cost_X_va)^.5*(cost_X_va/(PL*(1+tax_pl))) - 60*Y*(cost_Y/cost_Y_va)^.75*(cost_Y_va/PL) + 100 ⟂ PL
            mc_PK, -60*X*(cost_X/cost_X_va)^.5*(cost_X_va/(PK*(1+tax_pk))) - 40*Y*(cost_Y/cost_Y_va)^.75*(cost_Y_va/PK) + 100 ⟂ PK
        end)

        # Income Balance
        @constraint(mcp, ib, 
            CONS - 100*PL - 100*PK - 40*tax_pl*PL*X*(cost_X/cost_X_va)^.5*(cost_X_va/(PL*(1+tax_pl)))  - 60*tax_pk*PK*X*(cost_X/cost_X_va)^.5*(cost_X_va/(PK*(1+tax_pk))) ⟂ CONS
        )

fix(PW,1;force=true)

set_optimizer_attribute(mcp, "major_iteration_limit", 0)
optimize!(mcp)

object_dictionary(mcp)
Benchmarkresults_mcp = Dict(
    k => value.(v) for 
    (k, v) in object_dictionary(mcp) #if v isa AbstractArray{VariableRef}
)
mc_PX
value(mc_PX)
all_constraints(mcp; include_variable_in_set_constraints = true)

set_parameter_value(tax_pl,.5)
set_parameter_value(tax_pk,.5)

set_optimizer_attribute(mcp, "major_iteration_limit", 10000)
optimize!(mcp)
Counter1_results = Dict(
    k => value.(v) for 
    (k, v) in object_dictionary(mcp) #if v isa AbstractArray{VariableRef}
)# TXL = 0.2;

# TXK = 0.2;

### MPSGE.jl version
using MPSGE
M = MPSGEModel()

@parameters(M, begin
    tax_pl, 0, (description = "Tax on labor in sector X")
    tax_pk, 0, (description = "Tax on capital in sector X")
end)

@sectors(M, begin
    X
    Y
    W
end)

@commodities(M, begin
    PX
    PY
    PW
    PK
    PL
end)

@consumer(M, CONS)

@production(M, X, [t=0, s=.5, va=>s=1], begin
    @output(PX, 120, t)
    @input(PY, 20, s)
    @input(PL, 40, va, taxes=[Tax(CONS,tax_pl)])
    @input(PK, 60, va, taxes=[Tax(CONS,tax_pk)])
end)

@production(M, Y, [t=0, s=.75, va=>s=1], begin
    @output(PY, 120, t)
    @input(PX, 20, s)
    @input(PL, 60, va)
    @input(PK, 40, va)
end)

@production(M, W, [t=0, s=1], begin
    @output(PW, 200, t)
    @input(PX, 100, s)
    @input(PY, 100, s)
end)

@demand(M, CONS, begin
    @final_demand(PW, 200)
    @endowment(PL, 100)
    @endowment(PK, 100)
end)

fix(PW, 1)

solve!(M, cumulative_iteration_limit=0)
M.object_dict

M.jump_model
M
Benchmarkresults = generate_report(M)
MPSGE.market_clearance(PX)
value(MPSGE.market_clearance(M[:PX]))
value(PX)

set_value!(M[:tax_pl], 0.5)
set_value!(M[:tax_pk], 0.5)
solve!(M)
Counter1_results_M = generate_report(M)
