Live Programming and Computing with Python, R, Sage, Octave, Maxima, Singular, Gap, GP, HTML & Macaulay2

[Eli]

Choose any language from the dropdown menu, write some code in the editor (you can copy one of the example codes below and paste it into the editor) , and click the "Compute" button to run it (reload browser if you do not see the editor).

Example 0:

Code: [Select all] [Expand/Collapse]

1. import numpy as np
2. import matplotlib.pyplot as plt
3. from mpl_toolkits.mplot3d import Axes3D
4. from mpl_toolkits.mplot3d import proj3d
5.  
6. fig = plt.figure(figsize=(10,10)) #Define figure size
7. ax = fig.add_subplot(111, projection='3d')
8. plt.rcParams['legend.fontsize'] = 20
9.  
10. np.random.seed(4294967294) #Used the random seed for consistency
11.  
12. mu_vec_1 = np.array([0,0,0])
13. cov_mat_1 = np.array([[1,0,0],[0,1,0],[0,0,1]])
14. class_1_sample = np.random.multivariate_normal(mu_vec_1, cov_mat_1, 30).T
15. assert class_1_sample.shape == (3,30), "The matrix dimensions is not 3x30"
16.  
17. mu_vec_2 = np.array([1,1,1])
18. cov_mat_2 = np.array([[1,0,0],[0,1,0],[0,0,1]])
19. class_2_sample = np.random.multivariate_normal(mu_vec_2, cov_mat_2, 30).T
20. assert class_2_sample.shape == (3,30), "The matrix dimensions is not 3x30"
21.  
22. ax.plot(class_1_sample[0,:], class_1_sample[1,:], class_1_sample[2,:], 'o', markersize=10, color='green', alpha=1.0, label='Class 1')
23. ax.plot(class_2_sample[0,:], class_2_sample[1,:], class_2_sample[2,:], 'o', markersize=10, alpha=1.0, color='red', label='Class 2')
24.  
25. plt.title('Data Samples for Classes 1 & 2', y=1.04)
26. ax.legend(loc='upper right')
27. plt.savefig('Multivariate_distr.png', bbox_inches='tight')
28. plt.show()

Example one:

Code: [Select all] [Expand/Collapse]

1. import matplotlib.pyplot as plt
2. import seaborn as sns
3. # Load Dataset
4. df = sns.load_dataset('iris')
5. # Plot
6. plt.figure(figsize=(10,8), dpi= 80)
7. sns.pairplot(df, kind="reg", hue="species")
8. plt.show()

Code: [Select all] [Expand/Collapse]

1. import warnings
2. warnings.filterwarnings('ignore')
3.  
4. import seaborn as sns
5. import matplotlib.pyplot as plt
6. import pandas as pd
7. sns.set()
8. iris = sns.load_dataset("iris")
9. sns.pairplot(iris, hue='species', height=2.5)
10. textstr = 'Created at www.tssfl.com'
11. plt.gcf().text(0.36, 0.985, textstr, fontsize=14, color='green')
12. plt.show()

Code: [Select all] [Expand/Collapse]

1. import seaborn as sns
2. import matplotlib.pyplot as plt
3. import numpy as np
4. import pandas as pd
5.  
6. df = pd.DataFrame(np.random.randn(250, 4), columns=[f"Label%s{i+1}" % " " for i in range(4)])
7. df["Choice"] = np.random.choice([1., 0.], size=len(df))
8. sns.pairplot(df, vars=df.columns[:-1], hue="Choice")
9. plt.show()

Example Two:

Code: [Select all] [Expand/Collapse]

1. import matplotlib.pyplot as plt
2. names = ['group_a', 'group_b', 'group_c']
3. values = [1, 10, 100]
4.  
5. plt.figure(figsize=(9, 3))
6.  
7. plt.subplot(131)
8. plt.bar(names, values)
9. plt.subplot(132)
10. plt.scatter(names, values)
11. plt.subplot(133)
12. plt.plot(names, values)
13. plt.suptitle('Categorical Plotting')
14. plt.show()

Example Three:

Code: [Select all] [Expand/Collapse]

1. import matplotlib.pyplot as plt
2. import pandas as pd
3. import seaborn as sns
4. # Import Dataset
5.  
6. df = pd.read_csv("https://github.com/selva86/datasets/raw/master/mtcars.csv")
7.  
8. # Plot
9. plt.figure(figsize=(12,10), dpi= 80)
10. sns.heatmap(df.corr(), xticklabels=df.corr().columns, yticklabels=df.corr().columns, cmap='RdYlGn', center=0, annot=True)
11.  
12. textstr = 'Created at www.tssfl.com'
13. plt.gcf().text(0.64, 0.893, textstr, fontsize=14, color='darkblue')
14. # Decorations
15. plt.title('Correlogram of mtcars', fontsize=22, y=1)
16. plt.xticks(fontsize=12)
17. plt.yticks(fontsize=12)
18. plt.show()

Example Four:

Code: [Select all] [Expand/Collapse]

1. import matplotlib.pyplot as plt
2. import numpy as np
3. import pandas as pd
4. # Import Data
5. df = pd.read_csv("https://github.com/selva86/datasets/raw/master/mpg_ggplot2.csv")
6.  
7. # Prepare data
8. x_var = 'displ'
9. groupby_var = 'class'
10. df_agg = df.loc[:, [x_var, groupby_var]].groupby(groupby_var)
11. vals = [df[x_var].values.tolist() for i, df in df_agg]
12.  
13. # Draw
14. plt.figure(figsize=(16,9), dpi= 80)
15. colors = [plt.cm.Spectral(i/float(len(vals)-1)) for i in range(len(vals))]
16. n, bins, patches = plt.hist(vals, 30, stacked=True, density=False, color=colors[:len(vals)])
17.  
18. # Decoration
19. plt.legend({group:col for group, col in zip(np.unique(df[groupby_var]).tolist(), colors[:len(vals)])})
20. plt.title(f"Stacked Histogram of ${x_var}$ colored by ${groupby_var}$", fontsize=22)
21. plt.xlabel(x_var)
22. plt.ylabel("Frequency")
23. plt.ylim(0, 25)
24. #plt.xticks(ticks=bins[::3], labels=[round(b,1) for b in bins[::3]])
25. plt.show()

Output for Example One:

[Forbidden_Technology]

I have seen a number of new and very interesting stuffs, including Python integration, superb @Eli

[Eli]

More languages have been included: @Forbidden_Technology @RealityKing @Simulink and others.

[Eli]

Example One:

Run this R example code (From RDRR):

Code: [Select all] [Expand/Collapse]

1. require(graphics)
2.  
3. nus <- c(0:5, 10, 20)
4.  
5. x <- seq(0, 4, length.out = 501)
6. plot(x, x, ylim = c(0, 6), ylab = "", type = "n",
7.      main = "Bessel Functions  I_nu(x)")
8. for(nu in nus) lines(x, besselI(x, nu = nu), col = nu + 2)
9. legend(0, 6, legend = paste("nu=", nus), col = nus + 2, lwd = 1)
10.  
11. x <- seq(0, 40, length.out = 801); yl <- c(-.5, 1)
12. plot(x, x, ylim = yl, ylab = "", type = "n",
13.      main = "Bessel Functions  J_nu(x)")
14. abline(h=0, v=0, lty=3)
15. for(nu in nus) lines(x, besselJ(x, nu = nu), col = nu + 2)
16. legend("topright", legend = paste("nu=", nus), col = nus + 2, lwd = 1, bty="n")
17.  
18. ## Negative nu's --------------------------------------------------
19. xx <- 2:7
20. nu <- seq(-10, 9, length.out = 2001)
21. ## --- I() --- --- --- ---
22. matplot(nu, t(outer(xx, nu, besselI)), type = "l", ylim = c(-50, 200),
23.         main = expression(paste("Bessel ", I[nu](x), " for fixed ", x,
24.                                 ",  as ", f(nu))),
25.         xlab = expression(nu))
26. abline(v = 0, col = "light gray", lty = 3)
27. legend(5, 200, legend = paste("x=", xx), col=seq(xx), lty=1:5)
28.  
29. ## --- J() --- --- --- ---
30. bJ <- t(outer(xx, nu, besselJ))
31. matplot(nu, bJ, type = "l", ylim = c(-500, 200),
32.         xlab = quote(nu), ylab = quote(J[nu](x)),
33.         main = expression(paste("Bessel ", J[nu](x), " for fixed ", x)))
34. abline(v = 0, col = "light gray", lty = 3)
35. legend("topright", legend = paste("x=", xx), col=seq(xx), lty=1:5)
36.  
37. ## ZOOM into right part:
38. matplot(nu[nu > -2], bJ[nu > -2,], type = "l",
39.         xlab = quote(nu), ylab = quote(J[nu](x)),
40.         main = expression(paste("Bessel ", J[nu](x), " for fixed ", x)))
41. abline(h=0, v = 0, col = "gray60", lty = 3)
42. legend("topright", legend = paste("x=", xx), col=seq(xx), lty=1:5)
43.  
44.  
45. ##---------------  x --> 0  -----------------------------
46. x0 <- 2^seq(-16, 5, length.out=256)
47. plot(range(x0), c(1e-40, 1), log = "xy", xlab = "x", ylab = "", type = "n",
48.      main = "Bessel Functions  J_nu(x)  near 0\n log - log  scale") ; axis(2, at=1)
49. for(nu in sort(c(nus, nus+0.5)))
50.     lines(x0, besselJ(x0, nu = nu), col = nu + 2, lty= 1+ (nu%%1 > 0))
51. legend("right", legend = paste("nu=", paste(nus, nus+0.5, sep=", ")),
52.        col = nus + 2, lwd = 1, bty="n")
53.  
54. x0 <- 2^seq(-10, 8, length.out=256)
55. plot(range(x0), 10^c(-100, 80), log = "xy", xlab = "x", ylab = "", type = "n",
56.      main = "Bessel Functions  K_nu(x)  near 0\n log - log  scale") ; axis(2, at=1)
57. for(nu in sort(c(nus, nus+0.5)))
58.     lines(x0, besselK(x0, nu = nu), col = nu + 2, lty= 1+ (nu%%1 > 0))
59. legend("topright", legend = paste("nu=", paste(nus, nus + 0.5, sep = ", ")),
60.        col = nus + 2, lwd = 1, bty="n")
61.  
62. x <- x[x > 0]
63. plot(x, x, ylim = c(1e-18, 1e11), log = "y", ylab = "", type = "n",
64.      main = "Bessel Functions  K_nu(x)"); axis(2, at=1)
65. for(nu in nus) lines(x, besselK(x, nu = nu), col = nu + 2)
66. legend(0, 1e-5, legend=paste("nu=", nus), col = nus + 2, lwd = 1)
67.  
68. yl <- c(-1.6, .6)
69. plot(x, x, ylim = yl, ylab = "", type = "n",
70.      main = "Bessel Functions  Y_nu(x)")
71. for(nu in nus){
72.     xx <- x[x > .6*nu]
73.     lines(xx, besselY(xx, nu=nu), col = nu+2)
74. }
75. legend(25, -.5, legend = paste("nu=", nus), col = nus+2, lwd = 1)
76.  
77. ## negative nu in bessel_Y -- was bogus for a long time
78. curve(besselY(x, -0.1), 0, 10, ylim = c(-3,1), ylab = "")
79. for(nu in c(seq(-0.2, -2, by = -0.1)))
80.   curve(besselY(x, nu), add = TRUE)
81. title(expression(besselY(x, nu) * "   " *
82.                  {nu == list(-0.1, -0.2, ..., -2)}))
83.  

Example Two

Code: [Select all] [Expand/Collapse]

1. #Visualizing the diamonds dataset consisting of prices and quality information from about 54,000 diamonds,
2. #Included in the ggplot2 package, see https://www2.stat.duke.edu/courses/Fall15/sta112.01/post/hw/HW1.html
3. library(ggplot2)
4. library(dplyr)
5.  
6. ggplot(data = diamonds, mapping = aes(x = clarity)) + geom_bar(aes(fill = cut))

Example Three

Code: [Select all] [Expand/Collapse]

1. #Remove x-axis labels
2. library(ggplot2)
3. ggplot(data = diamonds, mapping = aes(x = color)) + geom_bar(aes(fill = cut))+
4.   theme(axis.title.x=element_blank(),
5.         axis.text.x=element_blank(),
6.         axis.ticks.x=element_blank())

[Eli]

Test this Discrete Fourier Transform MATLAB code (run using Octave):

Code: [Select all] [Expand/Collapse]

1. clear all, close all, clc
2. n =  1024;
3. w = exp(-i*2*pi/n);
4. [I,J] = meshgrid(1:n, 1:n);
5. DFT = w.^((I-1).*(J-1));
6. imagesc(real(DFT));
7. title("Discrete Fourier Transform Matrix")

The equivalent code in Python is:

Code: [Select all] [Expand/Collapse]

1. import numpy as np
2. import matplotlib.pyplot as plt
3. n =  1024
4. w = np.exp(-1j*2*np.pi/n)
5. j,k = np.meshgrid(np.arange(n), np.arange(n))
6. # Compute DFT
7. DFT = np.power(w,j*k)
8. # Real part
9. DFT = np.real(DFT)
10. plt.imshow(DFT)
11. plt.show()

[Eli]

Test these Python codes from here:

Code: [Select all] [Expand/Collapse]

1. import numpy as np
2. import scipy.special
3. import matplotlib.pyplot as plt
4.  
5. # create a figure window
6. fig = plt.figure(1, figsize=(9,8))
7.  
8. # create arrays for a few Bessel functions and plot them
9. x = np.linspace(0, 20, 256)
10. j0 = scipy.special.jn(0, x)
11. j1 = scipy.special.jn(1, x)
12. y0 = scipy.special.yn(0, x)
13. y1 = scipy.special.yn(1, x)
14. ax1 = fig.add_subplot(321)
15. ax1.plot(x,j0, x,j1, x,y0, x,y1)
16. ax1.axhline(color="grey", ls="--", zorder=-1)
17. ax1.set_ylim(-1,1)
18. ax1.text(0.5, 0.95,'Bessel', ha='center', va='top',
19.      transform = ax1.transAxes)
20.  
21. textstr = 'Created at www.tssfl.com'
22. plt.gcf().text(0.36, 0.92, textstr, fontsize=14, color='green')
23. # gamma function
24. x = np.linspace(-3.5, 6., 3601)
25. g = scipy.special.gamma(x)
26. g = np.ma.masked_outside(g, -100, 400)
27. ax2 = fig.add_subplot(322)
28. ax2.plot(x,g)
29. ax2.set_xlim(-3.5, 6)
30. ax2.axhline(color="grey", ls="--", zorder=-1)
31. ax2.axvline(color="grey", ls="--", zorder=-1)
32. ax2.set_ylim(-20, 100)
33. ax2.text(0.5, 0.95,'Gamma', ha='center', va='top',
34.      transform = ax2.transAxes)
35.  
36. # error function
37. x = np.linspace(0, 2.5, 256)
38. ef = scipy.special.erf(x)
39. ax3 = fig.add_subplot(323)
40. ax3.plot(x,ef)
41. ax3.set_ylim(0,1.1)
42. ax3.text(0.5, 0.95,'Error', ha='center', va='top',
43.      transform = ax3.transAxes)
44.  
45. # Airy function
46. x = np.linspace(-15, 4, 256)
47. ai, aip, bi, bip = scipy.special.airy(x)
48. ax4 = fig.add_subplot(324)
49. ax4.plot(x,ai, x,bi)
50. ax4.axhline(color="grey", ls="--", zorder=-1)
51. ax4.axvline(color="grey", ls="--", zorder=-1)
52. ax4.set_xlim(-15,4)
53. ax4.set_ylim(-0.5,0.6)
54. ax4.text(0.5, 0.95,'Airy', ha='center', va='top',
55.      transform = ax4.transAxes)
56.  
57. # Legendre polynomials
58. x = np.linspace(-1, 1, 256)
59. lp0 = np.polyval(scipy.special.legendre(0),x)
60. lp1 = np.polyval(scipy.special.legendre(1),x)
61. lp2 = np.polyval(scipy.special.legendre(2),x)
62. lp3 = np.polyval(scipy.special.legendre(3),x)
63. ax5 = fig.add_subplot(325)
64. ax5.plot(x,lp0, x,lp1, x,lp2, x,lp3)
65. ax5.axhline(color="grey", ls="--", zorder=-1)
66. ax5.axvline(color="grey", ls="--", zorder=-1)
67. ax5.set_ylim(-1,1.1)
68. ax5.text(0.5, 0.9,'Legendre', ha='center', va='top',
69.      transform = ax5.transAxes)
70.  
71. # Laguerre polynomials
72. x = np.linspace(-5, 8, 256)
73. lg0 = np.polyval(scipy.special.laguerre(0),x)
74. lg1 = np.polyval(scipy.special.laguerre(1),x)
75. lg2 = np.polyval(scipy.special.laguerre(2),x)
76. lg3 = np.polyval(scipy.special.laguerre(3),x)
77. ax6 = fig.add_subplot(326)
78. ax6.plot(x,lg0, x,lg1, x,lg2, x,lg3)
79. ax6.axhline(color="grey", ls="--", zorder=-1)
80. ax6.axvline(color="grey", ls="--", zorder=-1)
81. ax6.set_xlim(-5,8)
82. ax6.set_ylim(-5,10)
83. ax6.text(0.5, 0.9,'Laguerre', ha='center', va='top',
84.      transform = ax6.transAxes)
85.  
86. plt.show()

Code: [Select all] [Expand/Collapse]

1. import numpy as np
2. import matplotlib.pyplot as plt
3. from scipy.integrate import odeint
4.  
5. def f(y, t, params):
6.     theta, omega = y      # unpack current values of y
7.     Q, d, Omega = params  # unpack parameters
8.     derivs = [omega,      # list of dy/dt=f functions
9.              -omega/Q + np.sin(theta) + d*np.cos(Omega*t)]
10.     return derivs
11.  
12. # Parameters
13. Q = 2.0          # quality factor (inverse damping)
14. d = 1.5          # forcing amplitude
15. Omega = 0.65     # drive frequency
16.  
17. # Initial values
18. theta0 = 0.0     # initial angular displacement
19. omega0 = 0.0     # initial angular velocity
20.  
21. # Bundle parameters for ODE solver
22. params = [Q, d, Omega]
23.  
24. # Bundle initial conditions for ODE solver
25. y0 = [theta0, omega0]
26.  
27. # Make time array for solution
28. tStop = 200.
29. tInc = 0.05
30. t = np.arange(0., tStop, tInc)
31.  
32. # Call the ODE solver
33. psoln = odeint(f, y0, t, args=(params,))
34.  
35. # Plot results
36. fig = plt.figure(1, figsize=(8,8))
37.  
38. # Plot theta as a function of time
39. ax1 = fig.add_subplot(311)
40. ax1.plot(t, psoln[:,0])
41. ax1.set_xlabel('time')
42. ax1.set_ylabel('theta')
43.  
44. # Plot omega as a function of time
45. ax2 = fig.add_subplot(312)
46. ax2.plot(t, psoln[:,1])
47. ax2.set_xlabel('time')
48. ax2.set_ylabel('omega')
49.  
50. # Plot omega vs theta
51. ax3 = fig.add_subplot(313)
52. twopi = 2.0*np.pi
53. ax3.plot(psoln[:,0]%twopi, psoln[:,1], '.', ms=1)
54. ax3.set_xlabel('theta')
55. ax3.set_ylabel('omega')
56. ax3.set_xlim(0., twopi)
57.  
58. plt.tight_layout()
59. plt.show()

Code: [Select all] [Expand/Collapse]

1. import numpy as np
2. from scipy import fftpack
3. import matplotlib.pyplot as plt
4.  
5. width = 2.0
6. freq = 0.5
7.  
8. t = np.linspace(-10, 10, 101)   # linearly space time array
9. g = np.exp(-np.abs(t)/width) * np.sin(2.0*np.pi*freq*t)
10.  
11. dt = t[1]-t[0]       # increment between times in time array
12.  
13. G = fftpack.fft(g)   # FFT of g
14. f = fftpack.fftfreq(g.size, d=dt) # frequenies f[i] of g[i]
15. f = fftpack.fftshift(f)     # shift frequencies from min to max
16. G = fftpack.fftshift(G)     # shift G order to coorespond to f
17.  
18. fig = plt.figure(1, figsize=(8,6), frameon=False)
19. ax1 = fig.add_subplot(211)
20. ax1.plot(t, g)
21. ax1.set_xlabel('t')
22. ax1.set_ylabel('g(t)')
23.  
24. ax2 = fig.add_subplot(212)
25. ax2.plot(f, np.real(G), color='dodgerblue', label='real part')
26. ax2.plot(f, np.imag(G), color='coral', label='imaginary part')
27. ax2.legend()
28. ax2.set_xlabel('f')
29. ax2.set_ylabel('G(f)')
30.  
31. plt.show()

Code: [Select all] [Expand/Collapse]

1. import numpy as np
2. import scipy.optimize
3. import matplotlib.pyplot as plt
4.  
5. def tdl(x):
6.     y = 8./x
7.     return np.tan(x) - np.sqrt(y*y-1.0)
8.  
9. # Find true roots
10.  
11. rx1 = scipy.optimize.brentq(tdl, 0.5, 0.49*np.pi)
12. rx2 = scipy.optimize.brentq(tdl, 0.51*np.pi, 1.49*np.pi)
13. rx3 = scipy.optimize.brentq(tdl, 1.51*np.pi, 2.49*np.pi)
14. rx = np.array([rx1, rx2, rx3])
15. ry = np.zeros(3)
16. # print using a list comprehension
17. print('\nTrue roots:')
18. print('\n'.join('f({0:0.5f}) = {1:0.2e}'.format(x, tdl(x)) for x in rx))
19.  
20. # Find false roots
21.  
22. rx1f = scipy.optimize.brentq(tdl, 0.49*np.pi, 0.51*np.pi)
23. rx2f = scipy.optimize.brentq(tdl, 1.49*np.pi, 1.51*np.pi)
24. rx3f = scipy.optimize.brentq(tdl, 2.49*np.pi, 2.51*np.pi)
25. rxf = np.array([rx1f, rx2f, rx3f])
26. # print using a list comprehension
27. print('\nFalse roots:')
28. print('\n'.join('f({0:0.5f}) = {1:0.2e}'.format(x, tdl(x)) for x in rxf))
29.  
30. # Plot function and various roots
31.  
32. x = np.linspace(0.7, 8, 128)
33. y = tdl(x)
34. # Create masked array for plotting
35. ymask = np.ma.masked_where(np.abs(y)>20., y)
36.  
37. plt.figure(figsize=(6, 4))
38. plt.plot(x, ymask)
39. plt.axhline(color='black')
40. plt.axvline(x=np.pi/2., color="gray", linestyle='--', zorder=-1)
41. plt.axvline(x=3.*np.pi/2., color="gray", linestyle='--', zorder=-1)
42. plt.axvline(x=5.*np.pi/2., color="gray", linestyle='--', zorder=-1)
43. plt.xlabel(r'$x$')
44. plt.ylabel(r'$\tan x - \sqrt{(8/x)^2-1}$')
45. plt.ylim(-8, 8)
46.  
47. plt.plot(rx, ry, 'og', ms=5, label='true roots')
48.  
49. plt.plot(rxf, ry, 'xr', ms=5, label='false roots')
50. plt.legend(numpoints=1, fontsize='small', loc = 'upper right',
51.            bbox_to_anchor = (0.92, 0.97))
52. plt.tight_layout()
53. plt.show()

[Eli]

Subplots in Python and how to share axes x and y.

Test the codes below create subplots (Codes are taken from Matplotlib website).

Code: [Select all] [Expand/Collapse]

1. import numpy as np
2. import matplotlib.pyplot as plt
3. fig, (ax1, ax2) = plt.subplots(2, sharex='col',
4.                         gridspec_kw={'hspace':0})
5. fig.suptitle("Title")
6. data1 = np.array([1,2,3,4,5,6,7,8,9])
7. data2 = np.array([9,8,7,6,5,4,3,2,1])
8. ax1.plot(data1, color="red")
9. ax2.plot(data2, color="green")
10. plt.show()

Code: [Select all] [Expand/Collapse]

1. import matplotlib.pyplot as plt
2. import numpy as np
3. fig, axs = plt.subplots(2, 2, sharex='col', sharey='row',
4.                         gridspec_kw={'hspace': 0, 'wspace': 0})
5. (ax1, ax2), (ax3, ax4) = axs
6. fig.suptitle('Sharing x per column, y per row')
7. # Create example data
8. x = np.linspace(0, 2 * np.pi, 400)
9. y = np.sin(x ** 2)
10. ax1.plot(x, y)
11. ax2.plot(x, y**2, 'tab:orange')
12. ax3.plot(x + 1, -y, 'tab:green')
13. ax4.plot(x + 2, -y**2, 'tab:red')
14.  
15. for ax in axs.flat:
16.     ax.label_outer()
17. plt.show()

Code: [Select all] [Expand/Collapse]

1. import matplotlib.pyplot as plt
2. import numpy as np
3. fig, axs = plt.subplots(3, sharex=True, sharey=True, gridspec_kw={'hspace': 0})
4. fig.suptitle('Sharing both axes')
5. # Create example data
6. x = np.linspace(0, 2 * np.pi, 400)
7. y = np.sin(x ** 2)
8. axs[0].plot(x, y ** 2)
9. axs[1].plot(x, 0.3 * y, 'o')
10. axs[2].plot(x, y, '+')
11.  
12. # Hide x labels and tick labels for all but bottom plot.
13. for ax in axs:
14.     ax.label_outer()
15. plt.show()

Code: [Select all] [Expand/Collapse]

1. import matplotlib.pyplot as plt
2. import numpy as np
3. fig, axs = plt.subplots(3, sharex=True, sharey=True)
4. fig.suptitle('Sharing both axes')
5. # Create example data
6. x = np.linspace(0, 2 * np.pi, 400)
7. y = np.sin(x ** 2)
8. axs[0].plot(x, y ** 2)
9. axs[1].plot(x, 0.3 * y, 'o')
10. axs[2].plot(x, y, '+')
11. plt.show()

Code: [Select all] [Expand/Collapse]

1. import matplotlib.pyplot as plt
2. import numpy as np
3. fig, (ax1, ax2) = plt.subplots(2, sharex=True)
4. fig.suptitle('Aligning x-axis using sharex')
5. # Create example data
6. x = np.linspace(0, 2 * np.pi, 400)
7. y = np.sin(x ** 2)
8. ax1.plot(x, y)
9. ax2.plot(x + 1, -y)
10. plt.show()

Code: [Select all] [Expand/Collapse]

1. import matplotlib.pyplot as plt
2. import numpy as np
3. fig, axs = plt.subplots(2, 2)
4. # Create example data
5. x = np.linspace(0, 2 * np.pi, 400)
6. y = np.sin(x ** 2)
7. axs[0, 0].plot(x, y)
8. axs[0, 0].set_title('Axis [0,0]')
9. axs[0, 1].plot(x, y, 'tab:orange')
10. axs[0, 1].set_title('Axis [0,1]')
11. axs[1, 0].plot(x, -y, 'tab:green')
12. axs[1, 0].set_title('Axis [1,0]')
13. axs[1, 1].plot(x, -y, 'tab:red')
14. axs[1, 1].set_title('Axis [1,1]')
15.  
16. for ax in axs.flat:
17.     ax.set(xlabel='x-label', ylabel='y-label')
18.  
19. # Hide x labels and tick labels for top plots and y ticks for right plots.
20. for ax in axs.flat:
21.     ax.label_outer()
22. plt.show()

Code: [Select all] [Expand/Collapse]

1. import matplotlib.pyplot as plt
2. import numpy as np
3. fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2)
4. fig.suptitle('Sharing x per column, y per row')
5. # Create example data
6. x = np.linspace(0, 2 * np.pi, 400)
7. y = np.sin(x ** 2)
8. ax1.plot(x, y)
9. ax2.plot(x, y**2, 'tab:orange')
10. ax3.plot(x, -y, 'tab:green')
11. ax4.plot(x, -y**2, 'tab:red')
12.  
13. for ax in fig.get_axes():
14.     ax.label_outer()
15. plt.show()

[Eli]

Here is a nice Python code by Jake to animate the Lorenz system in 3D (Test it here):

Code: [Select all] [Expand/Collapse]

1. import numpy as np
2. from scipy import integrate
3.  
4. from matplotlib import pyplot as plt
5. from mpl_toolkits.mplot3d import Axes3D
6. from matplotlib.colors import cnames
7. from matplotlib import animation
8.  
9. N_trajectories = 20
10.  
11. def lorentz_deriv(params, t0, sigma=10., beta=8./3, rho=28.0):
12.     """Compute the time-derivative of a Lorentz system."""
13.     x, y, z = params
14.     return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]
15.  
16. # Choose random starting points, uniformly distributed from -15 to 15
17. np.random.seed(1)
18. x0 = -15 + 30 * np.random.random((N_trajectories, 3))
19.  
20. # Solve for the trajectories
21. t = np.linspace(0, 4, 1000)
22. x_t = np.asarray([integrate.odeint(lorentz_deriv, x0i, t)
23.                   for x0i in x0])
24.  
25. # Set up figure & 3D axis for animation
26. fig = plt.figure()
27. ax = fig.add_axes([0, 0, 1, 1], projection='3d')
28. ax.axis('off')
29.  
30. # choose a different color for each trajectory
31. colors = plt.cm.jet(np.linspace(0, 1, N_trajectories))
32.  
33. # set up lines and points
34. lines = sum([ax.plot([], [], [], '-', c=c)
35.              for c in colors], [])
36. pts = sum([ax.plot([], [], [], 'o', c=c)
37.            for c in colors], [])
38.  
39. # prepare the axes limits
40. ax.set_xlim((-25, 25))
41. ax.set_ylim((-35, 35))
42. ax.set_zlim((5, 55))
43.  
44. # set point-of-view: specified by (altitude degrees, azimuth degrees)
45. ax.view_init(30, 0)
46.  
47. # initialization function: plot the background of each frame
48. def init():
49.     for line, pt in zip(lines, pts):
50.         line.set_data([], [])
51.         line.set_3d_properties([])
52.  
53.         pt.set_data([], [])
54.         pt.set_3d_properties([])
55.     return lines + pts
56.  
57. # animation function.  This will be called sequentially with the frame number
58. def animate(i):
59.     # we'll step two time-steps per frame.  This leads to nice results.
60.     i = (2 * i) % x_t.shape[1]
61.  
62.     for line, pt, xi in zip(lines, pts, x_t):
63.         x, y, z = xi[:i].T
64.         line.set_data(x, y)
65.         line.set_3d_properties(z)
66.  
67.         pt.set_data(x[-1:], y[-1:])
68.         pt.set_3d_properties(z[-1:])
69.  
70.     ax.view_init(30, 0.3 * i)
71.     fig.canvas.draw()
72.     return lines + pts
73.  
74. # instantiate the animator.
75. #anim = animation.FuncAnimation(fig, animate, init_func=init, frames=500, interval=30, blit=True, save_count=1500)
76. anim = animation.FuncAnimation(fig, animate, frames=500, interval=30, blit=True, save_count=1500)
77.  
78. #Save as mp4. This requires mplayer or ffmpeg to be installed
79. anim.save('lorenz_attractor.mp4', fps=15, extra_args=['-vcodec', 'libx264'])
80.  
81. #Alternative saving options:
82. """
83. f = r"animation.mp4" #Change formats .avi, .mov,
84. writermp4 = animation.FFMpegWriter(fps=15)
85. anim.save(f, writer=writermp4)
86.  
87. #Or as .gif
88. f = r"animation.gif"
89. writergif = animation.PillowWriter(fps=15)
90. anim.save(f, writer=writergif)
91. """
92. plt.show()

You may also wish to check saving options here.


You can render the output in the 3D box with this code (see here), check the Matploltlib documentation too:

Code: [Select all] [Expand/Collapse]

1. import numpy as np
2. import matplotlib.pyplot as plt
3. from mpl_toolkits.mplot3d import Axes3D
4. from scipy.integrate import odeint as ODEint
5.  
6. N_trajectories = 20
7.  
8. def lorentz_deriv(params, t0, sigma=10., beta=8./3, rho=28.0):
9.     """Compute the time-derivative of a Lorentz system."""
10.     x, y, z = params
11.     return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]
12.  
13. x = np.linspace(0, 20, 1000)
14. y, z = 10.*np.cos(x), 10.*np.sin(x) # something simple
15.  
16. fig = plt.figure()
17. ax = fig.add_subplot(1,2,1,projection='3d')
18. ax.plot(x, y, z)
19.  
20. # now Lorentz
21. times = np.linspace(0, 4, 1000)
22.  
23. start_pts = 30. - 15.*np.random.random((20,3))  # 20 random xyz starting values
24.  
25. trajectories = []
26. for start_pt in start_pts:
27.     trajectory = ODEint(lorentz_deriv, start_pt, times)
28.     trajectories.append(trajectory)
29.  
30. ax = fig.add_subplot(1,2,2,projection='3d')
31. for trajectory in trajectories:
32.     x, y, z = trajectory.T  # transpose and unpack
33.     # x, y, z = zip(*trajectory)  # this also works!
34.     ax.plot(x, y, z)
35. plt.show()

Code: [Select all] [Expand/Collapse]

1. """
2. ========================================
3. Create 2D bar graphs in different planes
4. ========================================
5.  
6. Demonstrates making a 3D plot which has 2D bar graphs projected onto
7. planes y=0, y=1, etc.
8. """
9.  
10. from mpl_toolkits.mplot3d import Axes3D
11. import matplotlib.pyplot as plt
12. import numpy as np
13.  
14. fig = plt.figure()
15. ax = fig.add_subplot(111, projection='3d')
16. for c, z in zip(['r', 'g', 'b', 'y'], [30, 20, 10, 0]):
17.     xs = np.arange(20)
18.     ys = np.random.rand(20)
19.  
20.     # You can provide either a single color or an array. To demonstrate this,
21.     # the first bar of each set will be colored cyan.
22.     cs = [c] * len(xs)
23.     cs[0] = 'c'
24.     ax.bar(xs, ys, zs=z, zdir='y', color=cs, alpha=0.8)
25.    
26. ax.set_title("3D plotting")
27. ax.set_xlabel('X')
28. ax.set_ylabel('Y')
29. ax.set_zlabel('Z')
30.  
31. plt.show()

[Eli]

Test this machine learning code with scikit-learn from here:

Code: [Select all] [Expand/Collapse]

1. # Author: Tom Dupré la Tour
2. # License: BSD 3 clause
3.  
4. import numpy as np
5. import matplotlib.pyplot as plt
6.  
7. from sklearn.preprocessing import KBinsDiscretizer
8. from sklearn.datasets import make_blobs
9.  
10. print(__doc__)
11.  
12. strategies = ['uniform', 'quantile', 'kmeans']
13.  
14. n_samples = 200
15. centers_0 = np.array([[0, 0], [0, 5], [2, 4], [8, 8]])
16. centers_1 = np.array([[0, 0], [3, 1]])
17.  
18. # construct the datasets
19. random_state = 42
20. X_list = [
21.     np.random.RandomState(random_state).uniform(-3, 3, size=(n_samples, 2)),
22.     make_blobs(n_samples=[n_samples // 10, n_samples * 4 // 10,
23.                           n_samples // 10, n_samples * 4 // 10],
24.                cluster_std=0.5, centers=centers_0,
25.                random_state=random_state)[0],
26.     make_blobs(n_samples=[n_samples // 5, n_samples * 4 // 5],
27.                cluster_std=0.5, centers=centers_1,
28.                random_state=random_state)[0],
29. ]
30.  
31. figure = plt.figure(figsize=(14, 9))
32. i = 1
33. for ds_cnt, X in enumerate(X_list):
34.  
35.     ax = plt.subplot(len(X_list), len(strategies) + 1, i)
36.     ax.scatter(X[:, 0], X[:, 1], edgecolors='k')
37.     if ds_cnt == 0:
38.         ax.set_title("Input data", size=14)
39.  
40.     xx, yy = np.meshgrid(
41.         np.linspace(X[:, 0].min(), X[:, 0].max(), 300),
42.         np.linspace(X[:, 1].min(), X[:, 1].max(), 300))
43.     grid = np.c_[xx.ravel(), yy.ravel()]
44.  
45.     ax.set_xlim(xx.min(), xx.max())
46.     ax.set_ylim(yy.min(), yy.max())
47.     ax.set_xticks(())
48.     ax.set_yticks(())
49.  
50.     i += 1
51.     # transform the dataset with KBinsDiscretizer
52.     for strategy in strategies:
53.         enc = KBinsDiscretizer(n_bins=4, encode='ordinal', strategy=strategy)
54.         enc.fit(X)
55.         grid_encoded = enc.transform(grid)
56.  
57.         ax = plt.subplot(len(X_list), len(strategies) + 1, i)
58.  
59.         # horizontal stripes
60.         horizontal = grid_encoded[:, 0].reshape(xx.shape)
61.         ax.contourf(xx, yy, horizontal, alpha=.5)
62.         # vertical stripes
63.         vertical = grid_encoded[:, 1].reshape(xx.shape)
64.         ax.contourf(xx, yy, vertical, alpha=.5)
65.  
66.         ax.scatter(X[:, 0], X[:, 1], edgecolors='k')
67.         ax.set_xlim(xx.min(), xx.max())
68.         ax.set_ylim(yy.min(), yy.max())
69.         ax.set_xticks(())
70.         ax.set_yticks(())
71.         if ds_cnt == 0:
72.             ax.set_title("strategy='%s'" % (strategy, ), size=14)
73.  
74.         i += 1
75.  
76. plt.tight_layout()
77. plt.show()

[Eli]

Test Healpy with the following piece of codes:

Code: [Select all] [Expand/Collapse]

1. import numpy as np
2. import healpy as hp
3. import matplotlib.pyplot as plt
4. NSIDE = 32
5. m = np.arange(hp.nside2npix(NSIDE))
6. hp.mollview(m, title="Mollview image RING")
7. plt.show()

Code: [Select all] [Expand/Collapse]

1. import numpy as np
2. import healpy as hp
3. import matplotlib.pyplot as plt
4. NSIDE = 32
5. m = np.arange(hp.nside2npix(NSIDE))
6. hp.mollview(m, nest=True, title="Mollview image NESTED")
7. plt.show()

Code: [Select all] [Expand/Collapse]

1. import numpy as np
2. import healpy as hp
3. import matplotlib.pyplot as plt
4. NSIDE = 32
5. m = np.arange(hp.nside2npix(NSIDE))
6. hp.mollview(m, coord=['G','E'], title='Histogram equalized Ecliptic', unit='mK', norm='hist', min=-1,max=1, xsize=2000)
7. hp.graticule()
8. plt.show()