2

I have converted following code from IDL and I replaced Broyden function of IDL by broyden1 function of python, are both similar ?

import numpy as np
import spiceypy as spice
import pandas as pd
import scipy as sp
import math as mt
from scipy.optimize.nonlin import NoConvergence 
from scipy.interpolate import interp1d 
import matplotlib.pyplot as plt
from scipy import signal as sg

print('Start of Program.\n')

def nonrelbroyfunc(nonrelguessx): 
    global passdf,passc,passf,passbetaa,passvra,passvza,passdeltab,passvrb,passvzb,passzb,passra,passza,passgamma
    return [passdf * passc / passf - np.cos(passbetaa) * passvra - np.sin(passbetaa) * passvza + np.sin(passdeltab) * passvrb + np.cos(passdeltab) * passvzb + np.cos(passbetaa - nonrelguessx[1]) * passvra + np.sin(passbetaa - nonrelguessx[1]) * passvza - np.sin(passdeltab - nonrelguessx[0]) * passvrb - np.cos(passdeltab - nonrelguessx[0]) * passvzb, (-1) * passzb * np.sin(passdeltab - nonrelguessx[0]) - np.sqrt(passra**2 + passza**2) * np.sin(passbetaa - nonrelguessx[1] - passgamma)]

def relbroyfunc(relguessx): 
    global passdf,passc,passf,passbetaa,passvra,passvza,passdeltab,passvrb,passvzb,passzb,passra,passza,passgamma
    return [passdf / passf - (1 + passvrb * np.sin(passdeltab-relguessx[0]) / passc + passvzb * np.cos(passdeltab-relguessx[0]) / passc + 0.5 * passvbsqd / passc**2 - passub/passc**2) / (1 + passvra * np.cos(passbetaa-relguessx[1]) / passc + passvza * np.sin(passbetaa-relguessx[1]) / passc + 0.5 * passvasqd / passc**2 - passua/passc**2) + (1 + passvrb * np.sin(passdeltab) / passc + passvzb * np.cos(passdeltab) / passc + 0.5 * passvbsqd / passc**2 - passub/passc**2) / (1 + passvra * np.cos(passbetaa) / passc + passvza * np.sin(passbetaa) / passc + 0.5 * passvasqd / passc**2 - passua/passc**2),(-1) * passzb * np.sin(passdeltab - relguessx[0]) - np.sqrt(passra**2 + passza**2) * np.sin(passbetaa - relguessx[1] - passgamma)]

nconvergence =  1000
freq         =  8423*pow(10,6)

for iobs in range(0,nobs):
    print(iobs)
    passdf = deltaf[iobs]
    passf = freq
    passc = c_mks
    passbetaa = betaa
    passdeltab = deltab
    passvra = vra
    passvza = vza
    passvrb = vrb
    passvzb = vzb
    passzb = zb
    passra = ra
    passza = za
    passgamma = gamma

    nonrelguessx = [0,0] 
    try :
        nonrelresult = sp.optimize.broyden1( nonrelbroyfunc, nonrelguessx, maxiter = 5000, x_tol = 10**-20)
        converged = True
    except NoConvergence as e:
        nonrelresult = e.args[0]
        converged = False

    nonreldeltararr[iobs] = nonrelresult[0]
    nonrelbetararr[iobs] = nonrelresult[1]

    relguessx = [0,0] 
    try :
        relresult = sp.optimize.broyden1( relbroyfunc, relguessx, maxiter = 5000, x_tol = 10**-20)
        converged = True
    except NoConvergence as e:
        relresult = e.args[0]
        converged = False

    reldeltararr[iobs] = relresult[0]
    relbetararr[iobs] = relresult[1]

Values assigned to all the pass variables are derived from computation and those values are correct. I am trying to find nonreldeltararr, nonrelbetararr, reldeltararr and relbetararr. When applying broyden1 function , I am getting answer but not of the order I want. I think I am not properly applying broyden1 function or is there any another equivalent function in python to that of Broyden function in IDL. As there is no mention about Please, any help will be appreciated.

The documentation of broyden function in IDL

The documentation of broyden1 function in python

IDL code which I am implementing in python2.7 is as follows :

; Improve IDL's chances of using these functions correctly
; Make these two functions separate programs if problems are encountered
forward_function nonrelbroyfunc, relbroyfunc

; Common block used for solving a pair of equations for two critical angles
common passstuff, passdf, passf, passc, passbetaa, passdeltab, passvra, passvza, passvrb, passvzb, passzb, passra, passza, passgamma, passua, passub, passvasqd, passvbsqd

freq     = 8423.d6
nconvergence = 1000L

for iobs = 0, nobs-1 do begin                    
; Variables to pass to the Broyden root-finding functions
passdf = deltaf(iobs)
passf = freq
passc = c_mks
passbetaa = betaa
passdeltab = deltab
passvra = vra
passvza = vza
passvrb = vrb
passvzb = vzb
passzb = zb
passra = ra
passza = za
passgamma = gamma

nonrelguessx = [delta_r, beta_r] ; Either make use of prior non-relativistic result from Fjeldbo et al. (1971) numerical method
nonrelguessx = [0d,0d] ; Or ensure independent result
nonrelresult = broyden(nonrelguessx, 'nonrelbroyfunc', itmax=5000, /double, tolx=1d-20)
nonreldeltararr(iobs) = nonrelresult[0]
nonrelbetararr(iobs) = nonrelresult[1]

relguessx = [delta_r, beta_r] ; Either make use of prior non-relativistic result from Fjeldbo et al. (1971) numerical method
relguessx = [0d,0d] ; Or ensure independent result
relresult = broyden(relguessx, 'relbroyfunc', itmax=5000, /double, tolx=1d-20)
reldeltararr(iobs) = relresult[0]
relbetararr(iobs) = relresult[1]

endfor ; for iobs = 0, nobs-1 do begin

stop
end

FUNCTION nonrelbroyfunc, nonrelguessx 
common passstuff
RETURN, $
[passdf * passc / passf $ 
- cos(passbetaa) * passvra $
- sin(passbetaa) * passvza $
+ sin(passdeltab) * passvrb $
+ cos(passdeltab) * passvzb $
+ cos(passbetaa - nonrelguessx[1]) * passvra $
+ sin(passbetaa - nonrelguessx[1]) * passvza $
- sin(passdeltab - nonrelguessx[0]) * passvrb $
- cos(passdeltab - nonrelguessx[0]) * passvzb, $
(-1d) * passzb * sin(passdeltab - nonrelguessx[0]) - sqrt(passra^2 + passza^2) * sin(passbetaa - nonrelguessx[1] - passgamma)]
END 

FUNCTION relbroyfunc, relguessx 
common passstuff
RETURN, $
[passdf / passf - $
(1d + passvrb * sin(passdeltab-relguessx[0]) / passc + passvzb * cos(passdeltab-relguessx[0]) / passc + $
0.5d * passvbsqd / passc^2 - passub/passc^2) / $
(1d + passvra * cos(passbetaa-relguessx[1]) / passc + passvza * sin(passbetaa-relguessx[1]) / passc + $
0.5d * passvasqd / passc^2 - passua/passc^2) + $
(1d + passvrb * sin(passdeltab) / passc + passvzb * cos(passdeltab) / passc + $
0.5d * passvbsqd / passc^2 - passub/passc^2) / $
(1d + passvra * cos(passbetaa) / passc + passvza * sin(passbetaa) / passc + $
0.5d * passvasqd / passc^2 - passua/passc^2), $
(-1d) * passzb * sin(passdeltab - relguessx[0]) - sqrt(passra^2 + passza^2) * sin(passbetaa - relguessx[1] - passgamma)]
END 
Jajal
  • 21
  • 2

1 Answers1

1

There is an example in the IDL documentation page you linked to, have you tried to implement it in python? The bit of code below returns exactly the same values as in the example. This means that they are at least similar.

import numpy as np
import scipy.optimize

def fcn(x):
    return np.array([3.0 * x[0] - np.cos(x[1]*x[2]) - 0.5,
                    x[0]**2 - 81.0*(x[1] + 0.1)**2 + np.sin(x[2]) + 1.06,
                    np.exp(-x[0]*x[1]) + 20.0 * x[2] + (10.0*np.pi - 3.0)/3.0])

scipy.optimize.broyden1(fcn, np.array([1.0, 1.0, 1.0]))
# returns: array([  5.00000011e-01,  -2.34224116e-10,  -5.23598775e-01])

Edit: And since your bits of code are non-functional because incomplete we can't test them. But I see that relguessx is not defined in the python code.

ejb
  • 145
  • 8
  • Yes, I did tried it. Oh I see, its my mistake in writing down here but its correct in my original code and still it isn't working. – Jajal Mar 27 '18 at 07:06