/***********************************************
  Filename: example-05.c
  Purpose:
              Test GSL library.
              http://www.gnu.org/software/gsl/
  On Window XP & cygwin:
  On Ubuntu
        Compile: gcc -Wall -I/usr/local/include -c example-05.c
               Link: gcc -o example-05 -L/usr/local/lib example-05.o -lgsl -lgslcblas -lm

 On Mac OS X Lion:
        Compile: gcc -Wall -I/usr/local/include -c example-05.c
               Link: gcc -o example-05 -L/usr/local/lib example-05.o -lgsl -lgslcblas -lm

  Execute: ./example-05

  Result:
A is a 2 by 2 matrix
A =
[ [   3.00000000  4.00000000 ],
  [   1.00000000  2.00000000 ] ]
A =
[ [  -3.16227766 -4.42718872 ],
  [   0.16227766  0.63245553 ] ]
tau = [   1.94868330  0.00000000 ]
R =
[ [  -3.16227766 -4.42718872 ],
  [   0.00000000  0.63245553 ] ]
Q =
[ [   0.00000000  0.00000000 ],
  [   0.16227766  0.00000000 ] ]
***************************************************/

#include <stdio.h>
#include <string.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

void printVector(char *msg, gsl_vector *x, const char *fmt) {
   if (msg != NULL && strlen(msg) > 0)
       printf("%s", msg);
    size_t n = x->size;
    size_t i;
    printf("[ ");
    printf (fmt, gsl_vector_get(x, 0));
    for (i = 1; i < n; i++) {
        printf (fmt, gsl_vector_get(x, i));
    }
    printf(" ]\n");
}

void printMatrix(char *msg, gsl_matrix *a, char *fmt) {
    if (msg != NULL && strlen(msg) > 0)
        printf("%s\n", msg);
    size_t m = a->size1, n = a->size2;
    size_t i, j;
    printf("[ [ ");
    printf (fmt, gsl_matrix_get(a, 0, 0));
    for (j = 1; j < n; j++) {
        printf (fmt, gsl_matrix_get(a, 0, j));
    }
    printf(" ]");
    if (m > 1)
        printf(",\n");
    for (i = 1; i < m; i++) {
        printf("  [ ");
        printf (fmt, gsl_matrix_get(a, i, 0));
        for (j = 1; j < n; j++) {
            printf (fmt, gsl_matrix_get(a, i, j));
    }
    printf(" ]");
    if (i < m -1)
        printf(",\n");
    }
    printf(" ]\n");
}

void product(gsl_matrix *c, const gsl_matrix *a, const gsl_matrix *b) {
    size_t m = a->size1, n = a->size2, p = b->size2;
    size_t i, j, k;
    double sum;
    for (i = 0; i < m; i++) {
        for (k = 0; k < p; k++) {
            sum = 0.0;
            for (j = 0; j < n; j++) {
                sum = sum + gsl_matrix_get(a, i, j)*gsl_matrix_get(b, j, k);
            }
            gsl_matrix_set (c, i, k, sum);
        }
    }
}

void productWIthRightVector(gsl_vector *c, const gsl_matrix *a, const gsl_vector *b) {
    size_t m = a->size1, n = b->size;
    size_t i, j;
    double sum;
    for (i = 0; i < m; i++) {
        sum = 0.0;
        for (j = 0; j < n; j++) {
            sum = sum + gsl_matrix_get(a, i, j)*gsl_vector_get(b, j);
        }
        gsl_vector_set(c, i, sum);
    }
}

void productWIthLeftVector(gsl_vector *c, const gsl_vector *a, const gsl_matrix *b) {
    size_t m = a->size, n = b->size2;
    size_t i, j;
    double sum;
    for (j = 0; j < n; j++) {
        sum = 0.0;
        for (i = 0; i < m; i++) {
            sum = sum + gsl_vector_get(a, j)*gsl_matrix_get(b, i, j);
        }
        gsl_vector_set(c, j, sum);
    }
}

void add(gsl_matrix *c, const gsl_matrix *a, const gsl_matrix *b) {
    size_t m = a->size1, n = a->size2;
    size_t i, j;
    double v;
    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++) {
            v = gsl_matrix_get (a, i, j) + gsl_matrix_get (b, i, j);
         gsl_matrix_set (c, i, j, v);
  }
 }
}

void subtract(gsl_matrix *c, const gsl_matrix *a, const gsl_matrix *b) {
    size_t m = a->size1, n = a->size2;
    size_t i, j;
    double v;
    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++) {
            v = gsl_matrix_get (a, i, j) - gsl_matrix_get (b, i, j);
            gsl_matrix_set (c, i, j, v);
        }
    }
}

void multiply(gsl_matrix *c, double x, const gsl_matrix *a) {
    size_t m = a->size1, n = a->size2;
    size_t i, j;
    double v;
    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++) {
            v = x*gsl_matrix_get (a, i, j);
         gsl_matrix_set (c, i, j, v);
  }
 }
}

void set(gsl_matrix *c, const gsl_matrix *a) {
    gsl_matrix_memcpy(c, a);
}

double trace(const gsl_matrix *a) {
 size_t m = a->size1;
 size_t i;
 double sum = 0;
    for (i = 0; i < m; i++) {
        sum = sum + gsl_matrix_get (a, i, i);
 }
 return sum;
}

void transpose(gsl_matrix *c, const gsl_matrix *a) {
    size_t m = a->size1, n = a->size2;
    size_t i, j;
    for (i = 0; i < n; i++) {
        for (j = 0; j < m; j++) {
         gsl_matrix_set (c, i, j, gsl_matrix_get (a, j, i));
  }
 }
}

void getQR(gsl_matrix *Q, gsl_matrix *R, const gsl_matrix *A, const gsl_vector *tau) {
   size_t m = A->size1;
   size_t n = A->size2;
   size_t i, j;
   for (i =0 ; i < m; i++) {
       for (j =0 ; j < n; j++) {
           if (i > j) {
               gsl_matrix_set(R, i, j, 0);
               gsl_matrix_set(Q, i, j, gsl_matrix_get(A, i, j));
           }
           else {
               gsl_matrix_set(R, i, j, gsl_matrix_get(A, i, j));
               gsl_matrix_set(Q, i, j, 0);
           }
       }
   }
}

int main(void)  {
    size_t n = 2;
    gsl_vector *tau = gsl_vector_alloc(n);
    // size_t i, j;
    gsl_matrix *A = gsl_matrix_alloc(n, n);
    gsl_matrix *Q = gsl_matrix_alloc(n, n);
    gsl_matrix *R = gsl_matrix_alloc(n, n);

    gsl_matrix_set (A, 0, 0, 3);
    gsl_matrix_set (A, 0, 1, 4);
    gsl_matrix_set (A, 1, 0, 1);
    gsl_matrix_set (A, 1, 1, 2);

    printf ("A is a 2 by 2 matrix\n");
    printMatrix("A = ", A, "%12.8f");

    gsl_linalg_QR_decomp(A, tau);

    printMatrix("A = ", A, "%12.8f");
    printVector("tau = ", tau, "%12.8f");

    getQR(Q, R, A, tau);
    printMatrix("R = ", R, "%12.8f");
    printMatrix("Q = ", Q, "%12.8f");

     gsl_matrix_free (Q);
     gsl_matrix_free (R);
     gsl_matrix_free (A);
     gsl_vector_free (tau);
     return 0;
}

"""
 Filename: testLU2.py

    Find a LU decomposition of a given square matrix by using
    the function 'xalglib.rmatrixlu' of the library ALGLIB
    (www.alglib.net), Sergey Bochkanov and Vladimir Bystritsky.

    And find also the inverse of a given invertible square matrix
    by using the function 'xalglib.rmatrixinverse' of the library ALGLIN.

 Liscense: GPLv2
 Author: Copyright (c) 2011/09/19 (Mon)  PKim  pkim01 (AT) paran (DOT) com
"""

import math
import xalglib

class RVector:
    def __init__(self, arr):
        a = []
        for i in range(len(arr)):
            a.append(arr[i])
        self.arr = a
    def __str__(self):
        return str(self.arr)
    def __getitem__(self, key):
        return self.arr[key]
    def __setitem__(self, key, value):
        self.arr[i] = value
    def __add__(self, other):
        m = self.size()
        a = []
        for i in range(m):
           a.append(self[i] + other[i])
        return RVector(a)
    def __neg__(self):
        m = self.size()
        a = []
        for i in range(m):
           a.append(-self[i])
        return RVector(a)
    def __sub__(self, other):
        m = self.size()
        a = []
        for i in range(m):
           a.append(self[i] - other[i])
        return RVector(a)
    def __rmul__(self, other):
        if type(other) == int or type(other) == long or type(other) == float:
            m = self.size()
            a = []
            for i in range(m):
               a.append(other*self[i])
            return RVector(a)
        else:
            return None
    def angle(self, other):
        return math.acos( self.dot(other) / (self.dot(self) * other.dot(other)) )
    def cross(self, other):
        a = []
        value = self[1]*other[2] - self[2]*other[1]
        a.append(value)
        value = self[2]*other[0] - self[0]*other[2]
        a.append(value)
        value = self[0]*other[1] - self[1]*other[0]
        a.append(value)
        return RVector(a)
    def norm(self):
        n = self.size()
        sum = 0
        for i in range(n):
            sum = sum + self[i]*self[i]
        return math.sqrt(sum)
    def dot(self, other):
        n = self.size()
        sum = 0
        for i in range(n):
            sum = sum + self[i]*other[i]
        return sum
    def size(self):
        return len(self.arr)
    def printVector(self, msg):
        if msg != None and len(msg) > 0:
            print msg,
        print self

            Test fraction calculations, wgich works on Python 2.6 or 2.7.

   Execute: python TestFraction.py
      Date: 2011. 10. 10
"""

from fractions import *
import numpy as np
import scipy as sp
from scipy import *

def evalFraction(x):
 return x.numerator/float(x.denominator)

def printMatrix(msg, m):
 if msg != None and len(msg) > 0:
  print msg
 print "%s" % m

def minorMatrix(m, x, y):
 nr, nc = m.shape
 t = []
 for i in range(0, nr - 1):
  tmp = []
  for j in range(0, nc - 1):
   if i < x:
    if j < y:
     u = m[i, j]
     tmp.append(u)
    else:
     u = m[i, j+1]
     tmp.append(u)
   else:
    if j < y:
     u = m[i+1, j]
     tmp.append(u)
    else:
     u = Fraction(m[i+1, j+1])
     tmp.append(u)
  t.append(tmp)
 t2 = matrix(t)
 return t2
 
def det(m):
 nr, nc = m.shape
 if nr != nc:
      print "Caught TypeError"
      return None
 if nr == 1:
  return m[0, 0]
 elif nr == 2:
          return m[0, 0]*m[1, 1] - m[0, 1]*m[1, 0]
 else:
  sum = 0
  sign = 1
  for j in range(0, nc):
   minor = minorMatrix(m, 0, j)
   sum = sum + sign*m[0, j]*det(minor)
   sign = -sign
  return sum
 return None

def inverse(m):
 nr, nc = m.shape
 if nr != nc:
      print "Caught TypeError"
      return None
 if nr == 1:
  return matrix([[1/m[0, 0]]])
 elif nr == 2:
          d = [1, 1] - m[0, 1]*m[1, 0]
          return matrix( [[m[1, 1]/d, m[0,1]/d], [m[1, 0]/d, m[1, 1]/d]] )
 else:
  d = det(m)
  sign = 1
  t = []
  for j in range(0, nc):
   tmp = []
   for i in range(0, nr):
       minor = minorMatrix(m, i, j)
       sign = (-1)**(i + j)
       v = sign*det(minor)
       tmp.append(v/d)
   t.append(tmp)
  t2 = matrix(t)
  return t2
 return None

def hilbertMatrix(n):
 if n <= 0:
  return None
 t = []
 for i in range(0, n):
  tmp = []
  for j in range(0, n):
      tmp.append(Fraction(1, i + j + 1))
  t.append(tmp)
 t2 = matrix(t)
 return t2

z = Fraction(1, 2)
w = Fraction("5/4")

A = mat([[Fraction(1,5), Fraction(3,5)], [Fraction(3,5), Fraction(4,5)]])
B = mat([[Fraction(4,3), Fraction(-1,3)], [Fraction(2,3), Fraction(7,3)]])
printMatrix("A = ", A)
printMatrix("B = ", B)
printMatrix("A + B = ", A + B)
printMatrix("A - B = ", A - B)
printMatrix("A * B = ", A * B)
printMatrix("A / B = ", A / B)
printMatrix("A**2 = ", A**2)
printMatrix("A**3 = ", A**3)
printMatrix("B**-1 = ", B**-1)
printMatrix("B*B**-1 = ", B*B**-1)
printMatrix("B.transpose() = ", B.transpose())

print "B[0, 0] = %s" % (B[0, 0])
print "B[0, 1] = %s" % (B[0, 1])
print "B[1, 0] = %s" % (B[1, 0])
print "B[1, 1] = %s" % (B[1, 1])

B[0][0] = Fraction(5, 1) 
B[0][0] = 50
printMatrix("Changed: B = ", B)
print "  ---> det(B) = %s" % det(B)
B[0,0] = 80
printMatrix("Changed again: B = ", B)
print "  ---> det(B) = %s" % det(B)

C = mat([ [Fraction(1,1), Fraction(1,2), Fraction(1,3)],
          [Fraction(1,2), Fraction(1,3), Fraction(1,4)], 
          [Fraction(1,3), Fraction(1,4), Fraction(1,5)] ])
printMatrix("C = ", C)
v = det(C)
print "  ---> det(C) = %s" % v
print "       evalFraction(det(C)) = %s" % evalFraction(v)

from scipy import linalg
print "Cehck again by using scipy.linalg(): "
print "       linalg.det(C) = %s" % linalg.det(C)

D = inverse(C)
print "inverse(C) = \n%s" % D
printMatrix("inverse(C)*C = ", D*C)
printMatrix("C*inverse(C) = ", C*D)
printMatrix("linalg.inv(C) = ", linalg.inv(C))
printMatrix("linalg.inv(C)*C = ", linalg.inv(C)*C)

H = hilbertMatrix(4)
printMatrix("H = ", H)
v = det(H)
print "  ---> det(H) = %s" % v
print "       evalFraction(det(H)) = %s" % evalFraction(v)
print "Cehck again by using scipy.linalg(): "
print "       linalg.det(H) = %s" % linalg.det(H)
G = inverse(H)
print "inverse(H) = \n%s" % G
printMatrix("inverse(H)*H = ", G*H)
printMatrix("H*inverse(H) = ", H*G)
printMatrix("linalg.inv(H) = ", linalg.inv(H))
printMatrix("linalg.inv(H)*H = ", linalg.inv(H)*H)

"""
Execute Result:
A =
[[1/5 3/5]
 [3/5 4/5]]
B =
[[4/3 -1/3]
 [2/3 7/3]]
A + B =
[[23/15 4/15]
 [19/15 47/15]]
A - B =
[[-17/15 14/15]
 [-1/15 -23/15]]
A * B =
[[2/3 4/3]
 [4/3 5/3]]
A / B =
[[3/20 -9/5]
 [9/10 12/35]]
A**2 =
[[2/5 3/5]
 [3/5 1]]
A**3 =
[[11/25 18/25]
 [18/25 29/25]]
B**-1 =
[[ 0.7  0.1]
 [-0.2  0.4]]
B*B**-1 =
[[1.0 0.0]
 [-5.55111512313e-17 1.0]]
B.transpose() =
[[4/3 2/3]
 [-1/3 7/3]]
B[0, 0] = 4/3
B[0, 1] = -1/3
B[1, 0] = 2/3
B[1, 1] = 7/3
Changed: B =
[[50 50]
 [2/3 7/3]]
  ---> det(B) = 250/3
Changed again: B =
[[80 50]
 [2/3 7/3]]
  ---> det(B) = 460/3
C =
[[1 1/2 1/3]
 [1/2 1/3 1/4]
 [1/3 1/4 1/5]]
  ---> det(C) = 1/2160
       evalFraction(det(C)) = 0.000462962962963
Cehck again by using scipy.linalg():
       linalg.det(C) = 0.000462962962963
inverse(C) =
[[9 -36 30]
 [-36 192 -180]
 [30 -180 180]]
inverse(C)*C =
[[1 0 0]
 [0 1 0]
 [0 0 1]]
C*inverse(C) =
[[1 0 0]
 [0 1 0]
 [0 0 1]]
linalg.inv(C) =
[[   9.  -36.   30.]
 [ -36.  192. -180.]
 [  30. -180.  180.]]
linalg.inv(C)*C =
[[1.0 0.0 0.0]
 [0.0 1.0 0.0]
 [0.0 0.0 1.0]]
H =
[[1 1/2 1/3 1/4]
 [1/2 1/3 1/4 1/5]
 [1/3 1/4 1/5 1/6]
 [1/4 1/5 1/6 1/7]]
  ---> det(H) = 1/6048000
       evalFraction(det(H)) = 1.65343915344e-07
Cehck again by using scipy.linalg():
       linalg.det(H) = 1.65343915344e-07
inverse(H) =
[[16 -120 240 -140]
 [-120 1200 -2700 1680]
 [240 -2700 6480 -4200]
 [-140 1680 -4200 2800]]
inverse(H)*H =
[[1 0 0 0]
 [0 1 0 0]
 [0 0 1 0]
 [0 0 0 1]]
H*inverse(H) =
[[1 0 0 0]
 [0 1 0 0]
 [0 0 1 0]
 [0 0 0 1]]
linalg.inv(H) =
[[   16.  -120.   240.  -140.]
 [ -120.  1200. -2700.  1680.]
 [  240. -2700.  6480. -4200.]
 [ -140.  1680. -4200.  2800.]]
linalg.inv(H)*H =
[[1.0 0.0 0.0 -3.5527136788e-15]
 [0.0 1.0 -5.68434188608e-14 2.84217094304e-14]
 [0.0 0.0 1.0 1.13686837722e-13]
 [0.0 0.0 5.68434188608e-14 1.0]]
"""

local matrix = require "lua_matrix"

function printMatrix(msg, a)
    if type(a) == "number" then
        print( msg .. a )
        return
    end
    if type(msg) == "string" and #(msg) > 0 then
        print( msg .. "[ ")
    end
    print( a:tostring() .. " ]" )
end

function matrix.getRow(a, k)
    local n = #a[1]
    return matrix.subm(a, k, 1, k, n)
end

function matrix.row(a, k)
    local n = #a[1]
    return matrix.subm(a, k, 1, k, n)
end

function matrix.col(a, k)
    local n = #a
    return matrix.subm(a, 1, k, n, k)
end

local A, B, C, D

-- instantiate of matrices
A = matrix:new({{1, 2, 3}, {4, 5, 6}})
printMatrix( "A = ", A )
B = matrix:new({{3, 0, 1}, {1, 2, -2}})
printMatrix( "B = ", B )
printMatrix( "B + B = ", B + B )
printMatrix( "B - B = ", B - B )
printMatrix( "B:transpose() * B = ", B:transpose() * B )
printMatrix( "3 * B = ", 3 * B )
printMatrix( "0.1 * B = ", 0.1 * B )
local x, y = matrix.size(B)
print( "matrix.size(B) = ".. x .. "x" .. y)

C = matrix.copy( 2*B )
printMatrix( "C = ", C )
C = (3*B):copy()
printMatrix( "C = ", C )

printMatrix( "2*A = ", 2*A )
printMatrix( "3*B = ", 3*B )
printMatrix( "2*A + 3*B = ", 2*A + 3*B )
printMatrix( "2*A - 3*B = ", 2*A - 3*B )
print( "B:getelement(2,3) = ".. B:getelement(2,3))
print( "More simply, B[2][3] = ".. B[2][3])
print( "5*A:getelement(1,1) - B:getelement(2,3) = ".. (5*A:getelement(1,1) - B:getelement(2,3)))
print( "More simply, Simply, *A[1][1] - B[2][3] = ".. (5*A[1][1] - B[2][3]))
printMatrix( "matrix.subm(A, 1, 1, 1, 3) = ", matrix.subm(A, 1, 1, 1, 3))
printMatrix( "matrix.subm(B, 2, 1, 2, 3) = ", matrix.subm(B, 2, 1, 2, 3))
printMatrix( "matrix.row(B, 2) = ", matrix.row(B, 2))
printMatrix( "3*matrix.getRow(B, 2) = ", 3*matrix.getRow(B, 2))
printMatrix( "More simply, 3*matrix.row(B, 2) = ", 3*matrix.row(B, 2))
printMatrix( "matrix.col(B, 2) = ", matrix.col(B, 2))
printMatrix( "B/2 = ", B/2 )
printMatrix( "1/2*B = ", 1/2*B )

--[[
Execute Result:
A = [
1       2       3
4       5       6 ]
B = [
3       0       1
1       2       -2 ]
B + B = [
6       0       2
2       4       -4 ]
B - B = [
0       0       0
0       0       0 ]
B:transpose() * B = [
10      2       1
2       4       -4
1       -4      5 ]
3 * B = [
9       0       3
3       6       -6 ]
0.1 * B = [
0.3     0       0.1
0.1     0.2     -0.2 ]
matrix.size(B) = 2x3
C = [
6       0       2
2       4       -4 ]
C = [
9       0       3
3       6       -6 ]
2*A = [
2       4       6
8       10      12 ]
3*B = [
9       0       3
3       6       -6 ]
2*A + 3*B = [
11      4       9
11      16      6 ]
2*A - 3*B = [
-7      4       3
5       4       18 ]
B:getelement(2,3) = -2
More simply, B[2][3] = -2
5*A:getelement(1,1) - B:getelement(2,3) = 7
More simply, Simply, *A[1][1] - B[2][3] = 7
matrix.subm(A, 1, 1, 1, 3) = [
1       2       3 ]
matrix.subm(B, 2, 1, 2, 3) = [
1       2       -2 ]
matrix.row(B, 2) = [
1       2       -2 ]
3*matrix.getRow(B, 2) = [
3       6       -6 ]
More simply, 3*matrix.row(B, 2) = [
3       6       -6 ]
matrix.col(B, 2) = [
0
2 ]
B/2 = [
1.5     0       0.5
0.5     1       -1 ]
1/2*B = [
1.5     0       0.5
0.5     1       -1 ]
--]]

set term win
set parametric
set isosamples 80, 30
set hidden
set key below

set title "Parametric Sphere"
set urange [0:2*pi]
set vrange [-pi/2:pi/2]
set ztics nomirror -2.0, 0.5, 2.0
set view 45,50
a = 3
splot a*cos(u)*cos(v), a*sin(u)*cos(v), a*sin(v)

from mpl_toolkits.mplot3d.axes3d import Axes3D
import matplotlib.pyplot as plt

# imports specific to the plots in this example
import numpy as np
from matplotlib import cm

# Twice as wide as it is tall.
fig = plt.figure(figsize=plt.figaspect(0.7))

#---- First subplot
ax = fig.add_subplot(1, 1, 1, projection='3d')
U = np.arange(0, 2*np.pi, 0.1)
V = np.arange(0, 2*np.pi, 0.1)
U, V = np.meshgrid(U, V)
a = 3
X = a*np.cos(U)*np.cos(V)
Y = a*np.sin(U)*np.cos(V)
Z = a*np.sin(V)
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet,
        linewidth=0, antialiased=False)
ax.set_zlim3d(-2.01, 2.01)

fig.colorbar(surf, shrink=0.5, aspect=10)

plt.show()

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FixedLocator, FormatStrFormatter
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure()
ax = fig.gca(projection='3d')
U = np.arange(-2*np.pi, 2*np.pi, 0.2)
ulen = len(U)
V = np.arange(-2*np.pi, 2*np.pi, 0.2)
vlen = len(V)
U, V = np.meshgrid(U, V)
a = 5
b = 2
X = (a + b*np.cos(V))*np.cos(U)
Y = (a + b*np.cos(V))*np.sin(U)
Z = np.sin(V)

colortuple = ('y', 'b')
colors = np.empty(X.shape, dtype=str)
for v in range(vlen):
    for u in range(ulen):
        colors[u, v] = colortuple[(u + v) % len(colortuple)]

surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, facecolors=colors,
        linewidth=0, antialiased=False)

ax.set_zlim3d(-1, 1)
ax.w_zaxis.set_major_locator(LinearLocator(6))
ax.w_zaxis.set_major_formatter(FormatStrFormatter('%.03f'))

plt.show()

from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
U = np.arange(-2*np.pi, 2*np.pi, 0.2)
ulen = len(U)
V = np.arange(-2*np.pi, 2*np.pi, 0.2)
vlen = len(V)
U, V = np.meshgrid(U, V)
a = 5
b = 2
X = (a + b*np.cos(V))*np.cos(U)
Y = (a + b*np.cos(V))*np.sin(U)
Z = np.sin(V)
ax.plot_wireframe(X, Y, Z, rstride=3, cstride=3)

plt.show()

from mpl_toolkits.mplot3d.axes3d import Axes3D
import matplotlib.pyplot as plt

# imports specific to the plots in this example
import numpy as np
from matplotlib import cm

# Twice as wide as it is tall.
fig = plt.figure(figsize=plt.figaspect(0.7))

#---- First subplot
ax = fig.add_subplot(111, projection='3d')
U = np.arange(0, 2*np.pi, 0.1)
V = np.arange(0, 2*np.pi, 0.1)
U, V = np.meshgrid(U, V)
a = 5
b = 2
X = (a + b*np.cos(V))*np.cos(U)
Y = (a + b*np.cos(V))*np.sin(U)
Z = b*np.sin(V)
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet,
        linewidth=0, antialiased=False)
ax.set_zlim3d(-2.01, 2.01)

fig.colorbar(surf, shrink=0.5, aspect=10)

plt.show()