complex numbers

47 views
Skip to first unread message

MJ

unread,
Nov 15, 2010, 3:39:20 AM11/15/10
to algopy
Hi!!

It's Maria again... :D

First of all I have to tell you that the module it's working perfectly
and I get my matrices and the corresponding derivatives without a
problem.

But...of course...as I told you I'm using this program to model
optical systems, and therefore one of the main variables (let's say a
certain x[1], for example, in algopy) I have to use is the index of
refraction of the different materials light passes through, which is a
complex number. To test if algopy would work with them I set this
small example and I got an error.

>>> import numpy
>>> from numpy import matrix
>>> import cmath
>>> from cmath import sin, cos, sinh, cosh
>>> import algopy
>>> from algopy import UTPM, zeros
>>>
>>> def f(x):
... M = zeros((2,2), dtype=x)
... M[0,0] = M[1,1] = 1
... M[0,1] = sin(x[0]) + cos(x[1])
... M[1,0] = (sin(x[0]) + cos(x[1]))/2
... return M
...
>>> x = UTPM(zeros((2,2,2)))
>>> x.data[0,:] = [2+2j, 0+1j]
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
ValueError: setting an array element with a sequence.

Any ideas on how can I introduce this type of numbers in the process??

Thanks in advance!!!!

Maria.

Sebastian Walter

unread,
Nov 16, 2010, 3:58:01 AM11/16/10
to alg...@googlegroups.com
Hello Maria,

On Mon, Nov 15, 2010 at 9:39 AM, MJ <mlj.o...@gmail.com> wrote:
> Hi!!
>
> It's Maria again... :D
>
> First of all I have to tell you that the module it's working perfectly
> and I get my matrices and the corresponding derivatives without a
> problem.

I'm glad to hear that.

>
> But...of course...as I told you I'm using this program to model
> optical systems, and therefore one of the main variables (let's say a
> certain x[1], for example, in algopy) I have to use is the index of
> refraction of the different materials light passes through, which is a
> complex number. To test if algopy would work with them I set this
> small example and I got an error.

complex numbers are not (yet) officially supported.
Since there are no tests that include complex numbers,
it is possible that there are bugs in algopy that either may result in an error,
or worse, may result in an incorrect computation.

In conclusion: You can give it a try, but you'll have to be extra careful.
If you have small test problems that I could use to verify the
correctness of the computation,
I'd be glad to incorporate them into algopy.

>
>>>> import numpy
>>>> from numpy import matrix
>>>> import cmath
>>>> from cmath import sin, cos, sinh, cosh
>>>> import algopy
>>>> from algopy import UTPM, zeros
>>>>
>>>> def f(x):
> ...     M = zeros((2,2), dtype=x)
> ...     M[0,0] = M[1,1] = 1
> ...     M[0,1] = sin(x[0]) + cos(x[1])
> ...     M[1,0] = (sin(x[0]) + cos(x[1]))/2
> ...     return M
> ...
>>>> x = UTPM(zeros((2,2,2)))
>>>> x.data[0,:] = [2+2j, 0+1j]
> Traceback (most recent call last):
>  File "<stdin>", line 1, in <module>
> ValueError: setting an array element with a sequence.
>
> Any ideas on how can I introduce this type of numbers in the process??

I had to make some small modifications of algopy to get your example to work.
They will only be available in the next release 0.3.0.
If you want to give it a try now, you can get the bleeding edge
version from https://github.com/b45ch1/algopy ,
for instance on Linux with `$git clone git://github.com/b45ch1/algopy.git`.

You need to initialize the UTPM instance with a complex zeros array.
Below a code that should work

-------------------- start code ----------------------------
import numpy


import algopy
from algopy import UTPM, zeros

def f(x):
M = zeros((2,2), dtype=x)


M[0,0] = M[1,1] = 1

M[0,1] = numpy.sin(x[0]) + numpy.cos(x[1])
M[1,0] = (numpy.sin(x[0]) + numpy.cos(x[1]))/2
return M

x = UTPM(zeros((2,2,2),dtype=complex))


x.data[0,:] = [2+2j, 0+1j]

y = f(x)

-------------------- start code ----------------------------


In case you want the Jacobian of the function f you can use

----------------------- start code ------------------------
import numpy


import algopy
from algopy import UTPM, zeros

def f(x):
M = zeros((2,2), dtype=x)


M[0,0] = M[1,1] = 1

M[0,1] = numpy.sin(x[0]) + numpy.cos(x[1])
M[1,0] = (numpy.sin(x[0]) + numpy.cos(x[1]))/2
return M


x = UTPM.init_jacobian([2+2j, 0+1j])
y = f(x)
print UTPM.extract_jacobian(y)

----------------------- end code ------------------------

and get as output

python test_experimental.py[[[ 0.00000000+0.j 0.00000000+0.j ]
[-1.56562584-3.29789484j 0.00000000-1.17520119j]]

[[-0.78281292-1.64894742j 0.00000000-0.5876006j ]
[ 0.00000000+0.j 0.00000000+0.j ]]]

which is hopefully correct...


regards,
Sebastian

>
> Thanks in advance!!!!
>
> Maria.

MJ

unread,
Nov 16, 2010, 9:24:26 AM11/16/10
to algopy
Hey!!!!

wow...that's amazing!! you are really helpful..I can't wait until the
release of the new version... :D

On friday I have to show a raw version of my code so I'll use the "old
way" to perform the derivatives as a test. I'll check this results and
I'll try to implement algopy later. Just for clarity, the output is
the jacobian?? so the second example you wrote right? some of the
values seem to be according to my calculations but I don't understand
the order...I'll check it out and let you know.

THANK YOU VERY MUCH!!!

Maria

On Nov 16, 9:58 am, Sebastian Walter <sebastian.wal...@gmail.com>
wrote:
> Hello Maria,

Sebastian Walter

unread,
Nov 17, 2010, 4:19:21 AM11/17/10
to alg...@googlegroups.com
On Tue, Nov 16, 2010 at 3:24 PM, MJ <mlj.o...@gmail.com> wrote:
> Hey!!!!
>
> wow...that's amazing!! you are really helpful..I can't wait until the
> release of the new version... :D
>
> On friday I have to show a raw version of my code so I'll use the "old
> way" to perform the derivatives as a test. I'll check this results and
> I'll try to implement algopy later. Just for clarity, the output is
> the jacobian??
>so the second example you wrote right? some of the
> values seem to be according to my calculations but I don't understand
> the order...I'll check it out and let you know.

your function maps from R^2 --> R^{2 \times 2} and therefore the Jacobian is a
R^{2 \times 2 \times 2} matrix.

You can also reformulate your function such that it is a mapping
R^2 --> R^{4}

----------- start code -----------


import numpy
import algopy
from algopy import UTPM, zeros

def f(x):
M = zeros((2,2), dtype=x)
M[0,0] = M[1,1] = 1
M[0,1] = numpy.sin(x[0]) + numpy.cos(x[1])
M[1,0] = (numpy.sin(x[0]) + numpy.cos(x[1]))/2

return M.reshape(4)


x = UTPM.init_jacobian([2+2j, 0+1j])
y = f(x)
print UTPM.extract_jacobian(y)

------------ end code -------------------

then you get the R^{4 \times 2} Jacobian matrix as output:

MJ

unread,
Nov 17, 2010, 10:18:47 AM11/17/10
to algopy
Ok, it seems to be working perfectly fine because I get the same
values.

I haven't downloaded the test version because I have a mac and
apparently to use the "git..." command I have to download and install
some program (don't laugh but I haven't heard about this git-thing
before :), you see why I need the school?? ). I'll do it later and
I'll let you know how it works with my mirrors.

Thank you very much!!!...I'm learning a lot!!!

Maria.

On Nov 17, 10:19 am, Sebastian Walter <sebastian.wal...@gmail.com>
wrote:
Reply all
Reply to author
Forward
0 new messages