# Change the following if necessary:
Path = "C:/GISIranData/" ; Number_per_dot = 75
# Initializing variables
Points = 0 ; Crops = 0 ; Outlets = 0 ; Adjus = 0; Toler = 0
TotArea = 0  ; TotRev = 0 ; IterNo = 0 
TotDemand = [] # Dimensions - Crops 
Id = [] ; Area = [] ; Frac = [] # Dimensions - Points
Trans = [] # Dimensions - Outlets
Suit = [] ; Capac = [] ; Prod = [] ; ZoneRev = [] ; Select = [] # - Points X Crops
Price = [] ; NewPrice = [] # - Outlets X Crops
Dist = [] # - Points X Outlets
Rev = [] # - Crops X Outlets
# ----------------------
# 'Read' values from a .csv table
def read(Name):    
    global Crops,Points,Outlets,Adjus,Toler,TotDemand,Trans
    global Price,NewPrice, Path
    pTemp = [] ; import csv
    i = -1; stg = Path + Name + ".csv"
    with open(stg) as csvDataFile:
        csvReader = csv.reader(csvDataFile)       
        for row in csvReader:
            i = i+1
            if i == 1: 
                if Name == "Parameters":
                    Crops = int(row[0]) ;   Points = int(row[1])
                    Outlets = int(row[2]) ; Toler = int(row[3])
                    Adjus = float(row[4])
                if Name == "TotalDemand":
                    TotDemand.append(int(row[0]))
                    TotDemand.append(int(row[1]))
            if Name == "Outlets":
                for row in csvReader:
                    Trans.append(int(row[3])) ; jj = 3 ; pTemp = []
                    for j in range(Crops):
                        jj = jj + 1 ;  pTemp.append(int(row[jj]))
                    Price.append(pTemp) ; jj = 3 + Crops ; pTemp = []
                    for j in range(Crops):
                        jj = jj + 1 ;  pTemp.append(int(row[jj]))
                    NewPrice.append(pTemp)
            if Name == "Suitabilities":
                for row in csvReader:
                    stemp = [] ; l = 2
                    for j in range(0,Crops):
                        l = l + 1 ; stemp.append(int(row[l]))
                    Suit.append(stemp)
#--------------
# 'Observe' some variables in a layer
def observe(Name):
    global Area,TotArea,Frac,Dist,Capac,Id
    if Name == "Pattern":
        layer = QgsProject.instance().mapLayersByName(Name)[0]
        TotArea = 0
        for i in range(Points):
            feat = layer.getFeature(i)
            TotArea = TotArea + int(feat['Area'])
        for i in range(Points):
            feat = layer.getFeature(i)
            idtemp = [feat['Id']]; Id.append(idtemp)
            Area = int(feat['Area'])
            tempno = float((Area/TotArea) * Crops)
            Frac.append(tempno) ; DistTemp = [] ; tempstg = "Hello"
            for j in range(Outlets):
                tempstg = 'Dist_' + str(j+1) ; dtemp = int(feat[tempstg])
                DistTemp.append(dtemp)                   
            Dist.append(DistTemp) ; CapacTemp = []
            for j in range(Crops):
                ctemp = int(Frac[i] * TotDemand[j])
                ctemp = int((ctemp * (1 + ((Suit[i][j]-5)* 0.2))))
                CapacTemp.append(ctemp)                
            Capac.append(CapacTemp)
#------------ 
# Calculate optimal crop's production in each zone
def Find_soln():
    global Rev,TotRev,Prod,TotDemand,Adjus,Price,NewPrice,Toler
    global Crops,Outlets,Capac,Trans,Select
    Prod =[]; Select = [] ; TotRev = 0 ; Return = []
    for i in range(0,Points):
        Stemp = []
        for j in range(0,Crops):
            Stemp.append(0)
        Select.append(Stemp)
    for i in range(0,Points):
        Rev = []
        for j in range(0,Crops):
            OTemp = []
            for k in range(0,Outlets):
                rtemp = int(Capac[i][j] * NewPrice[k][j])
                rtemp = rtemp - int(Capac[i][j] * Trans[k])
                OTemp.append(rtemp)
            Rev.append(OTemp)
        M = - 100000 ; jj = 0 
        for j in range(0,Crops):
            for k in range(Outlets):
                if Rev[j][k] > M:
                    M = Rev[j][k]
        ties = 0 ; Return = [] ; MM = int(M - ((Toler/100) * M))
        for j in range(0,Crops):
            tempReturn = 0
            for k in range(Outlets):
                if Rev[j][k] > MM:
                    ties = ties + 1 ; Select[i][j] = 1
                    tempReturn = Rev[j][k] ; break 
            Return.append(tempReturn)
        ProdTemp = [] 
        for j in range(0,Crops):
            ptemp = 0
            if Select[i][j] == 1:
                TotRev = TotRev + int((Return[j]/ties)) 
                ptemp = int(Capac[i][j]/ties)
            ProdTemp.append(ptemp)
        Prod.append(ProdTemp) 
# Find consequent price changes
    for j in range(0,Crops):
        TotProd = 0
        for i in range(0,Points):
            TotProd = TotProd + Prod[i][j]
        print("Commodity",j+1," Demand =",TotDemand[j]," Prodn. =",TotProd)
        percentdiff = int(((TotProd - TotDemand[j]) * 100) / TotDemand[j])
        Miss = percentdiff * Adjus
        for i in range(0,Points):
            change = int((Miss / 100) * Capac[i][j])
            Capac[i][j] = int(Capac[i][j] - change)
        for k in range(0,Outlets):
            change = int((Miss / 100) * Price[k][j])
            NewPrice[k][j] = int(NewPrice[k][j] - change)
#------------
# Insert a new field (Name)
def New_fields(Layer_name,Name,Iter):
    global Path
    #tempstr = "C:\GISData\\" 
    tempstr = Path + Layer_name + ".shp" 
    layer = QgsVectorLayer(tempstr,"park","ogr")
    features=layer.getFeatures()
    from PyQt5.QtCore import QVariant
    layer_provider=layer.dataProvider()
    for k in range(Crops):
        tempstg = Name + str(k+1) + Iter
        layer_provider.addAttributes([QgsField(tempstg,QVariant.Int)])
        layer.updateFields()
        with edit(layer):
            for feature in layer.getFeatures():
                layer.updateFeature(feature)
        layer.commitChanges()
#------------
# Change values in a field (Name) to those in the (Variable) array
def Change_fields_vals(Layer_name,Name,Variable,Iter):
    global Prod, Path ; tempstr = Path + Layer_name + ".shp" 
    layer = QgsVectorLayer(tempstr,"park","ogr")
    features=layer.getFeatures()
    from PyQt5.QtCore import QVariant
    layer_provider=layer.dataProvider()
    with edit(layer):
        i = -1
        for feature in layer.getFeatures():
            i = i + 1
            for k in range(Crops):
                ii = k + 1 ; b = Name + str(ii) + Iter
                tnber = int(Variable[i][k]) ; feature[b] = tnber
                layer.updateFeature(feature)
# ----------
# Generates a new layer in which the optimal production of each
# crop in each zone is shown
def Record_optimal(Layer_name,Iter):
    global Prod, Path
    New_fields(Layer_name,"Prodn_",Iter)
    Change_fields_vals(Layer_name,"Prodn_",Prod,Iter)
    IterF = Iterations + 1
    if Iter == "_"+ str(IterF):
        tempstr = Path + Layer_name + ".shp"
        iface.addVectorLayer(tempstr,"Optimal", "ogr")
# ----------
# Make dot maps
def Make_dot_maps():
    import random ; global Prod,Crops,Path,Number_per_dot
    src = Path + "Pattern.shp"
    tractLyr = QgsVectorLayer(src, 'Pattern', "ogr")
    for i in range(0,Crops):
        xx = i + 1 ; tempstr = "Optimal_" + str(xx)
        popLyr =  QgsVectorLayer('Point?crs=epsg:4326', tempstr ,  "memory") 
        features = tractLyr.getFeatures() 
        vpr = popLyr.dataProvider() 
        j = -1
        for feature in features:
            j = j + 1 ; dotFeatures = [] 
            density = Prod[j][i] / Number_per_dot 
            found = 0 ; dots = [] ; g = feature.geometry() 
            minx =  g.boundingBox().xMinimum() 
            miny =  g.boundingBox().yMinimum() 
            maxx =  g.boundingBox().xMaximum() 
            maxy =  g.boundingBox().yMaximum()
            while found < density: 
                x = random.uniform(minx,maxx) 
                y = random.uniform(miny,maxy) 
                pnt = QgsPointXY(x,y)
                if g.contains(pnt):
                    dots.append(pnt) ; found += 1 
                geom = QgsGeometry.fromMultiPointXY(dots)
                f = QgsFeature() ; f.setGeometry(geom) 
                dotFeatures.append(f)
            vpr.addFeatures(dotFeatures) ; popLyr.updateExtents() 
        QgsProject.instance().addMapLayers([popLyr]) 
# --------------------------------------------------------------
# Run script
Iterations = 25 ; read("Parameters") ; read("TotalDemand") 
read("Outlets") ; read("Suitabilities") ; observe("Pattern") 
Find_soln() ; print("Total Revenue =",TotRev) ; ii = 1 
for i in range(Iterations):
    IterNo = ii ;  ii = ii + 1 ; print("Iteration",IterNo,":")
    Find_soln() ; print("Total Revenue =",TotRev)
    It = "_" + str(ii)
    if i == Iterations - 1:
        Record_optimal("Pattern",It)
Make_dot_maps()
# --------------------------------------------------------------
