$title The Farmer's Problem - Stochastic (FARMSP,SEQ=73)

$onText

This model helps a farmer to decide how to allocate
his or her land. The yields are uncertain.

Birge, R, and Louveaux, F V, Introduction to Stochastic Programming.
Springer, 1997.

$offText

Set crop  / wheat, corn, sugarbeets /
    ch    header for data table /
           yield  yield in tons per acre
           cost   plantcost on dollars per acre
           pprice crop seed purchase price in dollars per ton
           minreq minimum requirements of crop in ton to feed cattle /
alias (c,crop);

Table cd(crop,ch) crop data
           yield   cost  pprice  minreq
wheat        2.5    150     238     200
corn         3      230     210     240
sugarbeets  20      260
;

Parameter
   yf     yield factor / 1 /
   land   available land in acres /500/;

$onText
* Not clear what to do, that's why we do EMP!!!
Stochastic Parameter/Random Variable
    yf     yield factor / 1 /;
Set  r     realizations / below, avg, above  /;

Table yfdistrib(r,*)
        value prob
below     0.8 [1/4]
avg       1   [1/2]
above     1.2 [1/4];
option yf%yfdistrib;
yf.stage = 2;
yf.stage(t) = ord(t)+1;
$offText

Set seq price curve segments / s1*s2 /;
Table pricecurve(crop,seq,*) dollars per ton
              price     ub
wheat.s1        170    inf
corn.s1         150    inf
sugarbeets.s1    36   6000
sugarbeets.s2    10    inf
;
set pcs(crop,seq) relevant segments; option pcs<pricecurve;

set errorPC(crop) price curve is not concave;
errorPC(c) = smin(pcs(c,seq), pricecurve(c,seq,'price')-pricecurve(c,seq+1,'price'))<0;
abort$card(errorPC) errorPC;

Variables
   x(c)     crop planted in acres of land
   w(c,seq) crops sold in segment of cost curve in tons
   y(c)     crops purchased in tons
   profit   objective variable in dollars;
Positive variables x,w,y;

Equations
  profitdef  objective function
  landuse    capacity
  bal(c)     crop balance;

profitdef..    profit =e= sum(pcs, w(pcs)*pricecurve(pcs,'price'))
                        - sum(c, cd(c,'cost')*x(c) + cd(c,'pprice')*y(c));

landuse..      sum(c, x(c)) =l= land;

bal(c)..       yf*cd(c,'yield')*x(c) + y(c) - sum(pcs(c,seq), w(pcs)) =g= cd(c,'minreq');

* No purchase of crops that don't have a purchase price
y.fx(c)$(cd(c,'pprice')=0) = 0;
w.up(pcs) = pricecurve(pcs,'ub');

model farm_emp /all/;

file emp / '%emp.info%' /; put emp '* problem %gams.i%'/;
$onPut
randvar yf discrete 0.25 0.8
                    0.50 1.0
                    0.25 1.2

stage 2 yf y w bal profit
$offPut
putclose emp;

Set s            scenarios / s1*s3 /;
Parameter
    srep(s,*)    scenario attributes / #s.prob 0 /
    s_yf(s)      yield factor realization by scenario
    s_profit(s)  profit by scenario
    s_w(s,c,seq) crops sold in segment of cost curve in tons by scenario
    s_y(s,c)     crops purchased in tons by scenario;

Set dict / s     .scenario.''
           ''    .opt.     srep
           yf    .randvar. s_yf
           profit.level.   s_profit
           w     .level.   s_w
           y     .level.   s_y /;

$set EMPSOLVER %gams.emp%
$if x%EMPSOLVER% == x $set EMPSOLVER de

solve farm_emp using emp maximizing profit scenario dict;

display srep;

$if not set DEDICT $exit

$sTitle Generate dictionary for model generated by DE

set nodes / n1*n4 /;
option emp=de; farm_emp.optfile = 1;
file fdeopt / de.opt /; putclose fdeopt 'deDict farm_dict.txt' / 'subsolver convert';

solve farm_emp using emp maximizing profit scenario dict;
execute.checkErrorLevel 'rm -rf scratchdir && mkdir scratchdir && gams gams.gms lo=%gams.lo% scrdir=scratchdir';

Set extraUel / '', evobj, defevobj, defev,
                   varind, defvarind, defvartheta,
                   cvar, defcvar, cvardev, defcvardev, cvartarget, defcvartarget,
                   cvarobjdev, defcvarobjdev,
                   varcvarobj, varcvarind, defvarcvarind, defvarcvartheta /;
Equation profitdef_de(nodes), landuse_de(nodes), bal_de(c,nodes), xequ(*,*,*);
Variable x_de(c,nodes), w_de(c,seq,nodes), y_de(c,nodes), profit_de(nodes), yf_var_de(nodes), xvar(*,*,*);

embeddedCode Python:
from gams.core.gdx import *
import os
gams.debug = 1
def write_record(gdx, label, value, vals):
    vals[0] = float(value)
    rc = gdxDataWriteStr(gdx, [label], vals)
    assert rc == 1, f"gdxDataWriteStr (model) {label} failed"

def write_vesyms(gdx, xve_name, VE, XVE, uels, node_start, syms, keys, vals):
    last_sym_id = -1
    vals[0] = float(1)
    for ve in VE:
        if not int(ve[0]) == last_sym_id:
            if not last_sym_id == -1:
                rc = gdxDataWriteDone(gdx)
                assert rc == 1, f"gdxDataWriteDone {syms[last_sym_id]} failed"
            last_sym_id = int(ve[0])
            rc = gdxDataWriteRawStart(gdx, syms[last_sym_id]+'_de', "", len(ve)-1, GMS_DT_PAR , 0)
            assert rc == 1, f"gdxDataWriteRawStart {syms[last_sym_id]} failed"
        for i,k in enumerate(ve[1:]):
            keys[i] = int(k)+1
        keys[len(ve)-2] += node_start # node uel always last
        rc = gdxDataWriteRaw(gdx, keys, vals)
        assert rc == 1, f"gdxDataWriteRaw {syms[last_sym_id]} failed"
    rc = gdxDataWriteDone(gdx)
    assert rc == 1, f"gdxDataWriteDone {syms[last_sym_id]} failed"
    rc = gdxDataWriteStrStart(gdx, xve_name, "", 3, GMS_DT_PAR , 0)
    assert rc == 1, f"gdxDataWriteStrStart ({xve_name}) failed"
    for xve in XVE:
        l = [xve[0],'','']
        if len(xve)>1: l[1] = uels[int(xve[1])]
        if len(xve)>2: l[2] = uels[int(xve[2])]
        rc = gdxDataWriteStr(gdx, l, vals)
        assert rc == 1, f"gdxDataWriteStr ({xve_name}) failed"
    rc = gdxDataWriteDone(gdx)
    assert rc == 1, f"gdxDataWriteDone ({xve_name}) failed"

with open('farm_dict.txt') as f:
    nuel, nsym, xuels = [int(s) for s in f.readline().split()]
    uels = [ f.readline().strip() for i in range(nuel) ]
    syms = [ f.readline().strip() for i in range(nsym) ]
    equs = []
    vars = []
    xequs = []
    xvars = []
    ttlblk = 0
    while True:
        line = f.readline().strip()
        if not line: break
        if line[0:2] == 'EX':
            e = line.split()[1:]
            ttlblk += 3
            xequs.append(e)
        elif line[0] == 'E':
            e = line.strip().split()[1:]
            ttlblk += len(e)-1
            equs.append(e)
        elif line[0:2] == 'VX':
            v = line.strip().split()[1:]
            ttlblk += 3
            xvars.append(v)
        elif line[0] == 'V': # could be V or VR
            v = line.strip().split()[1:]
            ttlblk += len(v)-1
            vars.append(v)
        else:
            raise Exception(f'unknown line code of line >{line}<')

    gdxHandle = new_gdxHandle_tp()
    vals = doubleArray(5)
    keys = intArray(20)
    node_start = uels.index('n1')
    rc, msg = gdxCreate(gdxHandle, GMS_SSSIZE)
    if not rc:
        raise Exception(msg)
    rc = gdxOpenWrite(gdxHandle, os.path.join("scratchdir","gamsdict.%gams.scrExt%"), "DE dictionary")[0]
    assert rc == 1, "gdxOpenWrite failed"
    rc = gdxUELRegisterMapStart(gdxHandle)
    assert rc == 1, "gdxUELRegisterMapStart failed"
    for unr, u in enumerate(uels):
       rc = gdxUELRegisterMap(gdxHandle, unr+1, u)
       assert rc == 1, "gdxUELRegisterMap failed"
    rc = gdxUELRegisterDone(gdxHandle)
    assert rc == 1, "gdxUELRegisterDone failed"
    # Basic counts
    rc = gdxDataWriteStrStart(gdxHandle, "model", "basic counts", 1, GMS_DT_PAR , 0)
    assert rc == 1, "gdxDataWriteStrStart (model) failed"
    write_record(gdxHandle, 'ttlblk', ttlblk, vals)
    write_record(gdxHandle, 'mincolcnt', len(vars)+len(xvars), vals)
    write_record(gdxHandle, 'minrowcnt', len(equs)+len(xequs), vals)
    rc = gdxDataWriteDone(gdxHandle)
    assert rc == 1, "gdxDataWriteDone (model) failed"
    # Quoted uels, not necessary for gmscr (soution + dict -> GDX point file)
    rc = gdxDataWriteStrStart(gdxHandle, "display", "uels and their quotes", 1, GMS_DT_PAR , 0)
    assert rc == 1, "gdxDataWriteStrStart (display) failed"
    rc = gdxDataWriteDone(gdxHandle)
    assert rc == 1, "gdxDataWriteDone (display) failed"
    # Now write rows
    write_vesyms(gdxHandle, 'XEQU', equs, xequs, uels, node_start, syms, keys, vals)
    # Now write columns
    write_vesyms(gdxHandle, 'XVAR', vars, xvars, uels, node_start, syms, keys, vals)
    rc = gdxClose(gdxHandle)
    assert rc == 0, "gdxClose failed"
    gdxFree(gdxHandle)
endEmbeddedCode
* The tool gmscr takes gamssolu.dat (solution in vector format) and the dictionary (gamsdict.dat) and creates a GDX point file gmsgrid.gdx
$if     %system.filesys% == UNIX execute.checkErrorLevel 'gmscr_ux.out scratchdir%system.dirsep%gamscntr.%gams.scrExt%';
$if not %system.filesys% == UNIX execute.checkErrorLevel 'gmscr_nx.exe scratchdir%system.dirsep%gamscntr.%gams.scrExt%';
execute_loadpoint 'scratchdir%system.dirsep%gmsgrid.gdx',
    profitdef_de, landuse_de, bal_de, xequ, x_de, w_de, y_de, profit_de, yf_var_de, xvar;
display profitdef_de.l, landuse_de.m, bal_de.m, xequ.l;
display x_de.l, w_de.l, y_de.l, profit_de.l, yf_var_de.l, xvar.l;