Syntax conversion between Python and STATA

[dekoki]

Dear colleagues,

Im converting a code from Python to STATA that is used for coordinates transformation between two datums. I'm stuck at the point when have to apply the iterative procedure. Cannot find the way to write the loop, so that it includes the initial value and then swap the pointers.

Has any of you got any experience with that and could share some wisdom with me?

Here is this bit of the code:

#Lat is obtained by an iterative proceedure:
lat = arctan2(z_2,(p*(1-e2_2))) #Initial value
latold = 2*pi
while abs(lat - latold)>10**-16:
lat, latold = latold, lat
nu_2 = a_2/sqrt(1-e2_2*sin(latold)**2)
lat = arctan2(z_2+e2_2*nu_2*sin(latold), p)

Would appreciate any suggestions greatly,
Many thanks in advance!
De

[Dave Airey]

The only line of concern is the one with assignment to a target list "x,y = y,x", which flips values between two variables. Stata doesn't do this, so do it differently. Otherwise, Stata has while loops and math functions and operators.

[dekoki]

Ah, that was my worry.. Well, is there a way then I could inculde a little STATA programme within the loop, so that it would change the pointers and then proceed with iteration?

[Dave Airey]

One way is to use a temporary variable in your .do file as the third variable:

. do "/var/folders/86/_7btv8rj39q5pc1mc13_39_40000gn/T//SD03651.000000"

. clear

. set obs 10
obs was 0, now 10

. generate x = runiform()

. generate y = runiform() + 10

. list

     +---------------------+
     |        x          y |
     |---------------------|
  1. | .6039829   10.77213 |
  2. | .5235884   10.52387 |
  3. | .7738169   10.81054 |
  4. | .9424971   10.75908 |
  5. | .9735008   10.99192 |
     |---------------------|
  6. | .7693751   10.92293 |
  7. | .9032785   10.79384 |
  8. | .0548887   10.77065 |
  9. | .0163171   10.78728 |
 10. | .4746731   10.60474 |
     +---------------------+

. tempvar swap

. generate `swap' = x

. replace x = y
(10 real changes made)

. replace y = `swap'
(10 real changes made)

. list

     +--------------------------------+
     |        x          y   __000001 |
     |--------------------------------|
  1. | 10.77213   .6039829   .6039829 |
  2. | 10.52387   .5235884   .5235884 |
  3. | 10.81054   .7738169   .7738169 |
  4. | 10.75908   .9424971   .9424971 |
  5. | 10.99192   .9735008   .9735008 |
     |--------------------------------|
  6. | 10.92293   .7693751   .7693751 |
  7. | 10.79384   .9032785   .9032785 |
  8. | 10.77065   .0548887   .0548887 |
  9. | 10.78728   .0163171   .0163171 |
 10. | 10.60474   .4746731   .4746731 |
     +--------------------------------+

. 
end of do-file

. list

     +---------------------+
     |        x          y |
     |---------------------|
  1. | 10.77213   .6039829 |
  2. | 10.52387   .5235884 |
  3. | 10.81054   .7738169 |
  4. | 10.75908   .9424971 |
  5. | 10.99192   .9735008 |
     |---------------------|
  6. | 10.92293   .7693751 |
  7. | 10.79384   .9032785 |
  8. | 10.77065   .0548887 |
  9. | 10.78728   .0163171 |
 10. | 10.60474   .4746731 |
     +---------------------+

[Dave Airey]

Oh! There is swapval on SSC too. Do "ssc install swapval" in Stata. From looking at that, renaming and using a tempvar is probably better (efficient) than what I did or Nick would not have done it that way.

[dekoki]

Thanks, that helps!
I got to the point when im finally able to run it, but it seems like an infinite loop. It goes to the moment when I get a msg that there is no room for more variables.
Do you know what am I doing wrong here?

Here is the syntax:

*Lat is obtained by an iterative proceedure:
tempvar lat_2 latold
gen `lat_2' = atan2(Z_2,(P*(1-E2))) //Initial value latold
gen `latold' = 2*_pi

while abs(`lat'-`latold')>10^(-16){
tempvar lattemp
gen `lattemp' = `lat_2'
replace `lat_2' = `latold'
replace `latold' = `lattemp'
local NU = A/sqrt(1-E2*sin(`latold')^2)
local templat = atan2(Z_2+E2*NU*sin(`latold'), P)
}

[Robert Picard]

The only line of concern is the one with assignment to a target list "x,y = y,x", which flips values between two variables. Stata doesn't do this, so do it differently.

No need to copy variable content via a temporary variable, you can just do a group rename as in

rename (x y) (y x)

[Robert Picard]

Thanks, that helps!
I got to the point when im finally able to run it, but it seems like an infinite loop. It goes to the moment when I get a msg that there is no room for more variables.
Do you know what am I doing wrong here?

At every iteration, you are creating a new temporary variable and you do not drop it later. These variables accumulate until there is no more room to create new variables. As I explained in the previous post, this is not necessary as you can use a group rename instead.

You should also use doubles to store intermediate results. I don't think that using locals is right either. Try something like:

gen double lat = atan2(z_2, (p * (1 - e2_2)))
gen double latold = 2 * _pi

while abs(lat - latold) > 1e-16 {
rename (lat latold) (latold lat)
replace lat = atan2(z_2 + e2_2 * (a_2 / sqrt(1 - e2_2 * sin(latold)^2)) * sin(latold), p)
 }

[Dave Airey]

Robert, thanks! I forgot rename had been greatly enhanced, and the rename group help was worth reading just now. Thanks for jumping in and suggesting a better solution!

[Alan Riley (StataCorp)]

I got to the point when im finally able to run it, but it seems like an infinite loop. It goes to the moment when I get a msg that there is no room for more variables.
Do you know what am I doing wrong here?

Here is the syntax:

*Lat is obtained by an iterative proceedure:
tempvar lat_2 latold
gen `lat_2' = atan2(Z_2,(P*(1-E2))) //Initial value latold
gen `latold' = 2*_pi

while abs(`lat'-`latold')>10^(-16){
tempvar lattemp
gen `lattemp' = `lat_2'
replace `lat_2' = `latold'
replace `latold' = `lattemp'
local NU = A/sqrt(1-E2*sin(`latold')^2)
local templat = atan2(Z_2+E2*NU*sin(`latold'), P)
}

I haven't tried to run your code, but I see four potential problems.

1) This is actually two problems in one. I don't think the condition in your while loop is doing what you think it is. You coded

while abs(`lat'-`latold') > 10^(-16) {
   ...
}

The first problem is that I don't see `lat' anywhere else in your code. Do you mean `lat_2'?

The second problem is that what you coded is equivalent to

while abs(`lat'[1] - `latold'[1] > 10^-16) {
   ...
}

That is, your code is looking only at the values in the first observation of `lat' (which perhaps should be `lat_2') and `latold'. Are you trying to solve this for a vector of values or for a single value?


2) You should be using double variables for precision, e.g.

gen double `lattemp' = `lat_2'
...

3) You should be using scalars instead of locals for precision, e.g.

scalar NU = ...
scalar templat = ...

4) Speaking of that local templat line in your code, you don't refer to templat anywhere else in your code, so I don't think that's right. From your original code, I think you mean lat_2. If so, you need a replace command, not local.

[Robert Picard]

Alan is correct to point out the inconsistency of mixing scalar and variables in the code. I assumed that lat is a variable but I did not adjust the while condition. Here's one way to continue the while loop until all observations meet the threshold:

gen double lat = atan2(z_2, (p * (1 - e2_2)))
gen double latold = 2 * _pi

local more 1
while `more' {
rename (lat latold) (latold lat)
replace lat = atan2(z_2 + e2_2 * (a_2 / sqrt(1 - e2_2 * sin(latold)^2)) * sin(latold), p)
qui count if abs(lat - latold) > 1e-16
local more = r(N)
 }

[dekoki]

Dave, Robert, thanks very much for your input. Both solutions for swapping pointers work well!
Alan, you're right! It is supposed to be lat_2, not lat. I overlooked it somehow. The 'replace' option you and Robert proposed seems accurate, as well. Thanks for pointing this out!

gen double lat = atan2(z_2, (p * (1 - e2_2)))
gen double latold = 2 * _pi

while abs(lat - latold) > 1e-16 {
rename (lat latold) (latold lat)
replace lat = atan2(z_2 + e2_2 * (a_2 / sqrt(1 - e2_2 * sin(latold)^2)) * sin(latold), p)
}

This does magic, Robert. Thanks a lot!
I changed all my tempvars to double type and got expected results!
Many thanks, once again!

[Robert Picard]

I changed all my tempvars to double type and got expected results!

Glad to hear that you are getting the results you are expecting. Note that the code you cite in your last post uses an insufficient condition for the while loop. As Alan pointed out, this is equivalent to:

while abs(lat[1] - latold[1]) > 1e-16 {
...
 }

which means that only the values of lat and latold in the first observation are used to evaluate when to exit the loop. It may well be that when abs(lat[1] - latold[1]) > 1e-16 evaluates to false, all other observations are close enough to your expected results but you should check them all. Here's an example that illustrates the problem.

clear

set seed 123
set obs 10
gen double z_2 = uniform()
gen double e2_2 = uniform()
gen double p = uniform()
gen double a_2 = uniform()

* the following only tests the first observation

gen double lat = atan2(z_2, (p * (1 - e2_2)))
gen double latold = 2 * _pi

while abs(lat - latold) > 1e-16 {
rename (lat latold) (latold lat)
qui replace lat = atan2(z_2 + e2_2 * (a_2 / sqrt(1 - e2_2 * sin(latold)^2)) * sin(latold), p)
 }

gen delta = abs(lat - latold)
list lat latold delta


* redo with code that tests all observations

drop lat latold delta

gen double lat = atan2(z_2, (p * (1 - e2_2)))
gen double latold = 2 * _pi

local more 1
while `more' {
rename (lat latold) (latold lat)
qui replace lat = atan2(z_2 + e2_2 * (a_2 / sqrt(1 - e2_2 * sin(latold)^2)) * sin(latold), p)
qui count if abs(lat - latold) > 1e-16
local more = r(N)
 }

gen delta = abs(lat - latold)
list lat latold delta