Solving 3 non linear algebraic equations in python

[LarryBiodun]

How can I solve a non-linear algebraic equation in ArcGIS python over multiple rasters. Let me Rephrase. I want to solve the following 3 non linear equations ,    and for 46 8 day time steps.  I have 46 rasters each for an 8 day period for Β(σ) , and σ, where I need to take input values from per time step. Is/Io is a constant. The three unknowns I have are E, G and H. The solution should ideally produce 46 rasters each for the 3 unknowns. I don't know if this can be done in ArcGIS though. I know I can solve this equations easily in SciPy using the Fsolve function but I don't know how to get this working with rasters. Any code that could help will be greatly appreciated.

[LarryBiodun]

Hi Dan,

Thanks for your suggestion. But I honestly am struggling to fully understand it as I am still a newbie to python. I have included some sample data in the link below. If you have sometime. I'll appreciate if you can put this suggestion in a script as I really have no clue on how to write it. Thanks

https://drive.google.com/file/d/1dE32qKVm0V83d4pvmPUD0nWvDmrOaz3w/view?usp=sharing

[DanPatterson_Retired]

I looked.  The grids are floating point grids... there are millions of combinations, you will never pull out a lookup table nor will you be able to reclass the grids into something that combine will use.  Using a root finding solution just isn't going to work in this situation.  You might be advised to classify your input data into groupings of integer types, then come up with an representative value based on the combinations of the classes.... or... come up with an alternative approach, perhaps using different data.

[LarryBiodun]

Thanks for looking into it. Actually I found a possible way to do this, just that I technically dont know how to write the code. My brother has written something for me in C sharp but it is extremely slow because it is basically calling my python script and it is doing a lot of round tripping between the programs. It is in order of a day for a single set of files and I have 46 set of files.

The solution is as follows perhaps someone can help write a few lines in python to do this;

I converted the three input rasters to ASCII. Hence I have 4064 rows by 4876 columns. The program then calls each row and column at a time (i.e x,y), which essentially is a single value and computes fsolve for it. It then writes the answer for E, H and G from fsolve into 3 different .txt file for each x,y position. I  imagine this is doable in python as well, unfortunately, my brother doesnt know code in python and the non linear solvers in C Sharp have been problematic hence my conundrum. IF anyone can help write something along the lines of what I explained, it will be much appreciated.

import arcpy
import scipy
import numpy as np
def myFunction(z):
        E = z[0]
        G = z[1]
        H = z[2]
        F = np.empty((3))
        F[0] = (3.356/1.024)*3*H*pow(abs(H),-(1/6))-G
        F[1] = 300-E-H-G
        F[2] = 3.356*H-E
        return F

zGuess= np.array([1,1,1])
z = fsolve(myFunction,zGuess)

print z

[DanPatterson_Retired]

that is what I was saying would be ridiculously slow.

I suggested something similar a long time ago.  If you have 3 arrays/rasters, stacked one on top of the other, you pull your 3 values from each cell, giving you the 3 z's for your function.  then you would solve for that..

Give me 3 ascii files, one for each variable in ascii format, for the same area... how about 10 rows and 10 columns and I will see if I can simplify the solution for you.

I still think a reclassification to reduce the number of combinations is the way to go.  I just don't want to make up data

see this older thread  https://community.esri.com/message/732916-re-solving-3-non-linear-algebraic-equations-in-python?comm...

a0, b0 and c0 represent your three values for a location... Find like 'cells' that have that combination... chuck the tuple into the equation,throw the result into those cells.

Move on to the next unique triplet, repeat until you have run out of combinations.

Post data, it is doable, but I still think there must be a better way

[DanPatterson_Retired]

here is the principle... I have some work to do so you can fill in the technical gaps... I have kept it simple, I haven't included the determining the 'unique triplets', but If you are good with the idea and can test the principle on a small sample, I can fill in that bit

def dan_func(z):
    """feed in the input array, do the stuff and return"""
    d, e, f, ans = z.T
    ans = d + e + f
    z[:, 3] = ans
    return d, e, f, ans, z

a = np.array([[1, 4, 1],
              [1, 3, 1],
              [2, 4, 2]])

b = np.array([[1, 3, 1],
              [2, 4, 2],
              [1, 4, 1]])

c = np.array([[8, 7, 6],
              [6, 8, 8],
              [6, 8, 7]])

a0 = a.ravel()
b0 = b.ravel()
c0 = c.ravel()


z = np.zeros((9, 4))
z[:, 0] = a0
z[:, 1] = b0
z[:, 2] = c0

# do the work

d, e, f, ans, z = dan_func(z)
my_map = z[:, 3].reshape(a.shape)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Take 3 input rasters, unravel them, fill in an array that contains the number of cells wide by 1 more than the variables... ie a 3x3 raster with 4 variables gives you 'z ' which is 9x4

Call your function (dan_func is trademarked )

 Now this function simply adds the 3 values together and puts the result in the 4th column, then a bunch of stuff is returned, which you can print.  here is the 'z' array then it is reshaped into a map of 'dan stuff' defined by dan func.

z
Out[21]: 
array([[  1.,   1.,   8.,  10.],
       [  4.,   3.,   7.,  14.],
       [  1.,   1.,   6.,   8.],
       [  1.,   2.,   6.,   9.],
       [  3.,   4.,   8.,  15.],
       [  1.,   2.,   8.,  11.],
       [  2.,   1.,   6.,   9.],
       [  4.,   4.,   8.,  16.],
       [  2.,   1.,   7.,  10.]])

my_map
Out[22]: 
array([[ 10.,  14.,   8.],
       [  9.,  15.,  11.],
       [  9.,  16.,  10.]])‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

now tell your brother to give you a hand 

Good luck

[LarryBiodun]

Thanks Dan, I appreciate this. I am looking into this. I'll feed everyone back on my progress. Thanks all