MUI_IIMAI_1415_Introduccion

1096 days ago by trinidad1415

 

Sage es software libre y gratuito de matemáticas.

 

  • Libre: acceso al código fuente y licencia abierta   ⇨   podemos saber cómo funciona cada elemento, cambiarlo a nuestro gusto, distribuirlo(*)...

  • Gratuito (¿será bueno y bonito también?)

 

Érase una vez...

Sage fue creado en 2004-5 por William Stein, profesor de la U. de Washington. Su motivación: dejar de depender de programas de matemáticas caros y de código cerrado (¡un gran sinsentido en la investigación!). El objetivo: un programa con el que se pueda hacer matemáticas, equiparable a Magma, Mathematica, Matlab, Maple... en funcionalidad y eficiencia.

Dado todo lo que pueden hacer los programas comerciales, no era viable empezarlo desde cero. Por eso...

Sage es una interfaz entre muchos programas especializados en distintos aspectos de matemáticas, todos ellos libres y gratuitos. De hecho, Sage contiene todos esos programas en su instalación.

Como ocurre con algunos proyectos exitosos de código abierto, el desarrollo de Sage se basa en una comunidad muy fuerte.

 

Sage como interfaz

El lenguaje de programación elegido para ello es Python, un lenguaje moderno y relativamente sencillo de aprender.

Python es un lenguaje muy popular. Es usado por Google, Yahoo!, CERN, NASA, Industrial Light & Magic, etc., se incluye en casi todas las distribuciones de Linux...

⇨   Un alumno que aprende programación básica a través del uso de Sage está aprendiendo algo potencialmente muy útil.

Sage actúa de intermediario entre los distintos paquetes que incluye. Por otra parte, si estamos acostumbrados a usar algún programa que esté incluido en Sage podemos usar la sintaxis que ya conocemos. Y no es solo una interfaz: Sage está incorporando nuevas funciones constantemente, gracias a la comunidad de desarrolladores que tiene detrás.

 

 

Varias maneras de trabajar con Sage

  1. Conectado a un servidor donde tenga una cuenta, a través de un navegador (es decir, esto mismo)
    • La instalación local también es posible
  2. Por línea de comandos
  3. Usando una celda de Sage insertada en una página web, por ejemplo http://www.unex.es/eweb/sage/jornadamerida/ejemplocelda.html
    • Interacción con otros programas que "hablan" con páginas web, por ejemplo GeoGebra
  4. Si trabajamos con Python podemos importar Sage como librería
  5. Si usamos Latex podemos llamar a Sage para meter resultados de cálculos en nuestros documentos  ⇨  PDFs con hojas de ejercicios aleatorios, etc.

 

Las hojas de trabajo de Sage

Este documento es una hoja de trabajo. Consiste en celdas de cálculo (¡todavía no hemos visto ninguna!) y celdas de texto. Veamos una hoja que fue usada en una asignatura del CUM: http://sage.unex.es/home/pub/83

 

 

 

 

 

Vamos a dar una vuelta por mi cuenta de usuario...

 

 

 

 

¿ya de vuelta?

 

En resumen:

  • Hojas de trabajo individuales
  • Se pueden compartir con otros usuarios del mismo servidor
  • Se pueden publicar al mundo exterior
  • Se pueden exportar e importar

 

 

¡ Calculemos un poco !

 
1+1 
       
factorial(100) 
       
n(pi,1000) 
       
g = x^3-2*x^2+1 s = solve(g, x,solution_dict=True) s 
       
plot(g,(x,-1,2))+points([(k[x],0) for k in s],size=50,color='red') 
       

Cálculo diferencial e integral

f(x) = x * sin(2*x) 
       
plot(f(x),(x,-1,4),color='green',fill='axis', fillcolor='green')+plot(g(x)-g(0),(x,-1,4)) 
       

Álgebra lineal

A = matrix([[1,2],[4,5]]) A 
       
b = vector([1,-2]).column() b 
       
A^5 
       

¿$A$ es invertible?

A. 
       
A.inverse() 
       
A^(-1) 
       
v = A^(-1)*b v 
       
A*v == b 
       

Ecuaciones diferenciales

# orden dos con condicion inicial x = var('x') y = function('y', x) de = diff(y,x,2) - y == x desolve(de, y) 
       
f = desolve(de, y, [0,2,1]) print f plot(f,figsize=4) 
       
x = var('x') y = function('y', x) ff = [desolve(diff(y,x) + y - 1, y, [0,c]) for c in [-2,-1,0,1,2]] var('x y') plot_slope_field( -y + 1 , (x,0,2), (y,-2,2)) + plot(ff,(x,0,2),thickness=3) 
       

Programación básica

lista = range(3,30) lista 
       
[factor(i) for i in lista] 
       
[factor(2^i-1) for i in lista if is_prime(i)] 
       

Estadística

datos = [(80,174), (45,152), (63,160), (94,183), (24,102), (75,183), (56,148), (52,152), (61,166), (34,140), (21,98), (78,160)] datos 
       
X = [d[0] for d in datos] Y = [d[1] for d in datos] minX = min(X); maxX = -min([-t for t in X]) minY = min(Y); maxY = -min([-t for t in Y]) mX = mean(X); sX = std(X,bias=True) mY = mean(Y); sY = std(Y,bias=True) sXY = mean([(d[0]-mX)*(d[1]-mY) for d in datos]) (sXY/sX/sY).n() 
       
list_plot(datos,size=40,figsize=5) + line([(minX,mY),(maxX,mY)]) + line([(mX,minY),(mX,maxY)]) + plot(sXY/sX^2*(x-mX)+mY,(x,minX,maxX),color='red') 
       
alturas = [dd[1] for dd in datos] import matplotlib.pyplot as plt plt.figure(figsize=(6, 4)) plt.title(u"Histograma de alturas") plt.xlabel("metros") plt.ylabel("Frecuencia") plt.text(120, 1.6, r'$\mu=0.5$ (?)') plt.grid(True) plt.hist(alturas, facecolor='green') plt.savefig('Histogram.png') plt.close() 
       
# un test estadístico con Sage a,b = ttest([1,2,3,4,5],[1,2,3,3.5,5.121]) a; b 
       
ttest? 
       
# el mismo test llamando a R desde Sage en una línea r.t_test([1,2,3,4,5],[1,2,3,3.5,5.121]) 
       
%r # Toda la celda se ejecutará como si fuera R x <- c(1,2,3,4,5) y <- c(1,2,3,3.5,5.121) t.test(x,y) 
       
# Llamando a R desde Sage obtenemos un objeto con formato de R w = r.qnorm(0.025) w 
       
# Lo convertimos a Sage w.sage() 
       

(Lo bueno es que no tenemos por qué trabajar así, Sage hace esto "tras la cortina". Pero si alguna vez necesitamos hilar tan fino, podemos.)

 

Muestrario

A continuación veremos "de carrerilla" algunos ejemplos de lo que se puede hacer con Sage. En todos los casos podemos ver y cambiar el código o el contenido de la celda de texto correspondiente.

 

Fórmulas matemáticas

En esta frase metemos la fórmula $E=mc^2$ porque sí.

$$\sum_{\substack{
   0<i<m \\
   0<j<n
  }}
 P(i,j)$$

\[
 \lim_{x\to 0}{\frac{e^x-1}{2x}}
 \overset{\left[\frac{0}{0}\right]}{\underset{\mathrm{H}}{=}}
 \lim_{x\to 0}{\frac{e^x}{2}}={\frac{1}{2}}
\]



 

Insertar imágenes y vídeos

 

Antes hemos visto varias imágenes insertadas, también podemos insertar por ejemplo un vídeo de Youtube.

 

 

 

Interacciones

 

Con Sage es posible crear objetos interactivos para ilustrar conceptos, permitir la experimentación... en cada ejemplo que sigue encontrarás en el código (oculto hasta que se pulsa en "hide") el enlace de donde se ha copiado. Puedes ver estos y más en http://wiki.sagemath.org/interact/ y en http://interact.sagemath.org/top-rated-posts

g = graphs.HouseGraph() graph_editor(g); 
       
a = randint(0,5) b = randint(0,5) html('<h2>¿ Cuánto es '+str(a)+'+'+str(b)+' ?</h2>') @interact def _(y=(0..10),): if y == a+b: print '¿ '+str(a)+'+'+str(b)+' = '+str(y)+' ? BIEN' else: print '¿ '+str(a)+'+'+str(b)+' = '+str(y)+' ? MAL' 
       

Click to the left again to hide and once more to show the dynamic interactive window

#%hide html('<h2>Tangente a una curva en un punto</h2>') # http://wiki.sagemath.org/interact/calculus#A_simple_tangent_line_grapher @interact def tangent_line(f = input_box(default=sin(x)), xbegin = slider(0,10,1/10,0), xend = slider(0,10,1/10,10), x0 = slider(0, 1, 1/100, 1/2, label="i/d")): prange = [xbegin, xend] x0i = xbegin + x0*(xend-xbegin) var('x') df = diff(f) x0in = x0i.n() tanf = f(x=x0i) + df(x=x0i)*(x-x0i) tanfn = expand(f(x=x0in) + df(x=x0in)*(x-x0in)) fplot = plot(f, prange[0], prange[1]) print 'Tangent line is y = ' + tanf._repr_() print ' = ' + tanfn._repr_() tanplot = plot(tanf, prange[0], prange[1], rgbcolor = (1,0,0)) fmax = f.find_local_maximum(prange[0], prange[1])[0] fmin = f.find_local_minimum(prange[0], prange[1])[0] show(fplot + tanplot + point((x0i,f(x=x0i)),size=20), xmin = prange[0], xmax = prange[1], ymax = fmax, ymin = fmin, figsize=5) 
       

Click to the left again to hide and once more to show the dynamic interactive window

#%hide html('<h2>Integración aproximada con sumas de Riemann</h2>') # http://interact.sagemath.org/node/45 @interact def demo(g=sin(x),a=0,b=pi,n=slider(0,20,1,6),method=selector(['Left', 'Right', 'Midpoint', 'Trapezoid'], nrows=1, label="Method"),fillcurve=('Fill area under curve',False)): var('k') f(x)=g dx = (b-a)/n pts = [a+i*dx for i in range(0,n+1)] p=Graphics() if method=='Left': for i in range(0,n): h=f(pts[i]) p+=polygon([(pts[i],0),(pts[i],h),(pts[i+1],h),(pts[i+1],0)],fill=False) approx=dx*add([f(pts[i]) for i in range(0,n)]) if method=='Right': for i in range(0,n): h=f(pts[i+1]) p+=polygon([(pts[i],0),(pts[i],h),(pts[i+1],h),(pts[i+1],0)],fill=False) approx=dx*add([f(pts[i]) for i in range(1,n+1)]) if method=='Midpoint': midpts=[a+1/2*dx+i*dx for i in range(0,n+1)] for i in range(0,n): h=f(midpts[i]) p+=polygon([(pts[i],0),(pts[i],h),(pts[i+1],h),(pts[i+1],0)],fill=False) approx=dx*add([f(midpts[i]) for i in range(0,n)]) if method=='Trapezoid': for i in range(0,n): p+=polygon([(pts[i],0),(pts[i],f(pts[i])),(pts[i+1],f(pts[i+1])),(pts[i+1],0)],fill=False) approx=dx*add([(f(pts[i])+f(pts[i+1]))/2 for i in range(0,n)]) p+=plot(f(x),(x,a,b),color='red',fill=fillcurve) show(p) actual=integral(f(x),x,a,b).n() approx=approx.n() html(r'$\displaystyle\int_{%s}^{%s} f(x)\, dx = %s$' % (latex(a) , latex(b),actual)) html('<br/>') html(method +' Approximation = %s'% approx) html('Error = %s'% abs(actual-approx)) 
       

Click to the left again to hide and once more to show the dynamic interactive window

#%hide # http://wiki.sagemath.org/interact/linear_algebra#The_Gauss-Jordan_method_for_inverting_a_matrix html("<h2>Inversión de matrices por Gauss-Jordan</h2>") #Choose the size D of the square matrix: D = 3 example = [[1 if k==j else 0 for k in range(D)] for j in range(D)] example[0][-1] = 2 example[-1][0] = 3 @interact def _(M=input_grid(D,D, default = example, label='Matrix to invert', to_value=matrix), tt = text_control('Enter the bits of precision used' ' (only if you entered floating point numbers)'), precision = slider(5,100,5,20), auto_update=False): if det(M)==0: print 'Failure: Matrix is not invertible' return if M.base_ring() == RR: M = M.apply_map(RealField(precision)) N=M M=M.augment(identity_matrix(D)) print 'We construct the augmented matrix' show(M) for m in range(0,D-1): if M[m,m] == 0: lista = [(M[j,m],j) for j in range(m,D)] maxi, c = max(lista) M[c,:],M[m,:]=M[m,:],M[c,:] print 'We permute rows %d and %d'%(m+1,c+1) show(M) for n in range(m+1,D): a=M[m,m] if M[n,m]!=0: print "We add %s times row %d to row %d"%(-M[n,m]/a, m+1, n+1) M=M.with_added_multiple_of_row(n,m,-M[n,m]/a) show(M) for m in range(D-1,-1,-1): for n in range(m-1,-1,-1): a=M[m,m] if M[n,m]!=0: print "We add %s times row %d to the row %d"%(-M[n,m]/a, m+1, n+1) M=M.with_added_multiple_of_row(n,m,-M[n,m]/a) show(M) for m in range(0,D): if M[m,m]!=1: print 'We divide row %d by %s'%(m+1,M[m,m]) M = M.with_row_set_to_multiple_of_row(m,m,1/M[m,m]) show(M) M=M.submatrix(0,D,D) print 'We keep the right submatrix, which contains the inverse' html('$$M^{-1}=%s$$'%latex(M)) print 'We check it actually is the inverse' html('$$M^{-1} \cdot M=%s\cdot%s=%s$$'%(latex(M),latex(N),latex(M*N))) 
       

Click to the left again to hide and once more to show the dynamic interactive window

#%hide # http://wiki.sagemath.org/interact/misc#An_Interactive_Venn_Diagram html("<h2>Diagrama de Venn interactivo</h2>") def f(s, braces=True): t = ', '.join(sorted(list(s))) if braces: return '{' + t + '}' return t def g(s): return set(str(s).replace(',',' ').split()) @interact def _(X='1,2,3,a', Y='2,a,3,4,apple', Z='a,b,10,apple'): S = [g(X), g(Y), g(Z)] X,Y,Z = S XY = X & Y XZ = X & Z YZ = Y & Z XYZ = XY & Z # html('<center>') html(" $X \cap Y$ = %s"%f(XY)) html(" $X \cap Z$ = %s"%f(XZ)) html(" $Y \cap Z$ = %s"%f(YZ)) html("$X \cap Y \cap Z$ = %s"%f(XYZ)) # html('</center>') centers = [(cos(n*2*pi/3), sin(n*2*pi/3)) for n in [0,1,2]] scale = 1.7 clr = ['yellow', 'blue', 'green'] G = Graphics() for i in range(len(S)): G += circle(centers[i], scale, rgbcolor=clr[i], fill=True, alpha=0.3) for i in range(len(S)): G += circle(centers[i], scale, rgbcolor='black') # Plot what is in one but neither other for i in range(len(S)): Z = set(S[i]) for j in range(1,len(S)): Z = Z.difference(S[(i+j)%3]) G += text(f(Z,braces=False), (1.5*centers[i][0],1.7*centers[i][1]), rgbcolor='black') # Plot pairs of intersections for i in range(len(S)): Z = (set(S[i]) & S[(i+1)%3]) - set(XYZ) C = (1.3*cos(i*2*pi/3 + pi/3), 1.3*sin(i*2*pi/3 + pi/3)) G += text(f(Z,braces=False), C, rgbcolor='black') # Plot intersection of all three G += text(f(XYZ,braces=False), (0,0), rgbcolor='black') # Show it G.show(aspect_ratio=1, axes=False) 
       

Click to the left again to hide and once more to show the dynamic interactive window

#%hide # http://wiki.sagemath.org/interact/games#Zeros import random def init(): global B,Br,n,round,tm,t,v B = 40 Br = 15 n = 1 round = 0 tm = 0 t = walltime() init() html("<h2>¿Cuántos ceros hay?</h2>") @interact def zeros(a=selector(buttons=True, nrows=1, values=['Reset'] + [1..B], default=1)): if a == 'Reset': init() html("<center>") global B,Br,n,round,tm,t,v if a == n: if round > 0: html("<font size=+3 color='red'>RIGHT</font>") r = walltime() - t tm += r round += 1 t = walltime() while True: n2 = random.randrange(1,Br) if n2 != n: n = n2 break if Br < B: Br += 2 elif round > 0: html("<font size=+2 color='blue'>Wrong. Try again...</font>") html("</center>") html("<font size=+%s color='#333'>"%random.randrange(-2,5)) print ' '*random.randrange(20) + '0'*n html("</font>") if round > 0: html("<br><br><center>Score: %s rounds, Average time: %.1f seconds</center>"%( round, float(tm)/round)) 
       

Click to the left again to hide and once more to show the dynamic interactive window

#%hide # http://interact.sagemath.org/node/85 html("<h2>Método Diffie-Hellman de intercambio de claves (teoría de números)</h2>") @interact def diffie_hellman(bits=slider(8, 513, 4, 8, 'Number of bits', False), button=selector(["Show new example"],label='',buttons=True)): maxp = 2 ^ bits p = random_prime(maxp) k = GF(p) if bits > 100: g = k(2) else: g = k.multiplicative_generator() a = ZZ.random_element(10, maxp) b = ZZ.random_element(10, maxp) html(""" <style> .gamodp, .gbmodp { color:#000; padding:5px } .gamodp { background:#846FD8 } .gbmodp { background:#FFFC73 } .dhsame { color:#000; font-weight:bold } </style> <h2 style="color:#000;font-family:Arial, Helvetica, sans-serif">%s-Bit Diffie-Hellman Key Exchange</h2> <ol style="color:#000;font-family:Arial, Helvetica, sans-serif"> <li>Alice and Bob agree to use the prime number p = %s and base g = %s.</li> <li>Alice chooses the secret integer a = %s, then sends Bob (<span class="gamodp">g<sup>a</sup> mod p</span>):<br/>%s<sup>%s</sup> mod %s = <span class="gamodp">%s</span>.</li> <li>Bob chooses the secret integer b=%s, then sends Alice (<span class="gbmodp">g<sup>b</sup> mod p</span>):<br/>%s<sup>%s</sup> mod %s = <span class="gbmodp">%s</span>.</li> <li>Alice computes (<span class="gbmodp">g<sup>b</sup> mod p</span>)<sup>a</sup> mod p:<br/>%s<sup>%s</sup> mod %s = <span class="dhsame">%s</span>.</li> <li>Bob computes (<span class="gamodp">g<sup>a</sup> mod p</span>)<sup>b</sup> mod p:<br/>%s<sup>%s</sup> mod %s = <span class="dhsame">%s</span>.</li> </ol> """ % (bits, p, g, a, g, a, p, (g^a), b, g, b, p, (g^b), (g^b), a, p, (g^ b)^a, g^a, b, p, (g^a)^b)) 
       

Click to the left again to hide and once more to show the dynamic interactive window

#%hide html('<h2>Teorema del bocadillo (geometría plana)</h2>') # http://interact.sagemath.org/node/82 """Ham creates an interact in sage which lets the user input two non-intersecting polygons as ordered pairs and also specify an error bound. A line is generated to approximate the line which simultaneously bisects both areas of the polygons. This is displayed graphically with the equation of the line as well as the quality of the approximation of both bisections.""" #### Max and Min methods def find_ymax(ls1, ls2): """Given two lists of 2-tuples, ls1 and ls2, outputs the largest y-coordinate.""" Q=[] W=[] for i in range(len(ls1)): g=ls1[i] Q.append(g[1]) for i in range(len(ls2)): g=ls2[i] W.append(g[1]) ymax1 = max(Q) ymax2 = max(W) return max([ymax1, ymax2]) def find_ymin(ls1, ls2): """Given a list of 2-tuples, ls1 and ls2, outputs the smallest y-coordinate""" Q=[] W=[] for i in range(len(ls1)): g=ls1[i] Q.append(g[1]) for i in range(len(ls2)): g=ls2[i] W.append(g[1]) ymin1 = min(Q) ymin2 = min(W) return min(ymin1, ymin2) def find_xmax(ls1, ls2): """Given a list of 2-tuples, ls1 and ls2, outputs the largest x-coordinate""" Q=[] W=[] for i in range(len(ls1)): g=ls1[i] Q.append(g[0]) for i in range(len(ls2)): g=ls2[i] W.append(g[0]) xmax1 = max(Q) xmax2 = max(W) return max(xmax1, xmax2) def find_xmin(ls1, ls2): """Given a list of 2-tuples, ls1 and ls2, outputs the smallest x-coordinate""" Q=[] W=[] for i in range(len(ls1)): g=ls1[i] Q.append(g[0]) for i in range(len(ls2)): g=ls2[i] W.append(g[0]) xmin1 = min(Q) xmin2 = min(W) return min(xmin1, xmin2) def max_list(ls1, ls2): """Given two lists of 2-tuples, ls1 and ls2, outputs the length of the larger of the two lists.""" S=[] R=[] for i in range(0,len(ls1)): for z in ls1[i]: S.append(z) newS = S[0:len(S):2] for i in range(0,len(ls2)): for z in ls2[i]: R.append(z) newR = R[0:len(R):2] m=max(newR) n=max(newS) return max(m,n) def min_list(ls1, ls2): """Given two lists of tuples, ls1 and ls2, outputs the length of the smaller of the two lists.""" S=[] R=[] for i in range(0,len(ls1)): for z in ls1[i]: S.append(z) newS = S[0:len(S):2] for i in range(0,len(ls2)): for z in ls2[i]: R.append(z) newR = R[0:len(R):2] m=min(newR) n=min(newS) return min(m,n) #### Line Construction and Methods def draw_line(pt1, pt2): """Given two 2-tuple points, pt1 and pt2, outputs a line as a 3-tuple for ax+by=c.""" if pt2[0]-pt1[0] != 0: m = (pt2[1] - pt1[1])/(pt2[0]- pt1[0]) b = pt1[1] -m*(pt1[0]) newline = line() newline.a = -m newline.b = 1 newline.c = b return newline else: m = (pt2[0] - pt1[0])/(pt2[1]- pt1[1]) b = pt1[0] - m*(pt1[1]) newline = line() newline.a = 1 newline.b = -m newline.c = b return newline def give(self): """Converts a line to a list [a,b,c] for ax + by = x.""" ls = [self.a, self.b, self.c] return ls ####Polygon Class class pgon: def __init__(self, list_vertices = []): self.vertex_list = list_vertices self.area = self.PArea() self.xcent = self.cent_x() self.ycent = self.cent_y() self.centroid = (self.xcent, self.ycent) def PArea(self): """Computes area of polygon given a list of vertices.""" S=[] for i in range(0,len(self.vertex_list)): for z in self.vertex_list[i]: S.append(z) A=0 for i in range(0,len(S),2): if i < len(S)-2: A= A + S[i]*S[i+3] - S[i+2]*S[i+1] else: A= A + S[i]*S[1] - S[0]*S[i+1] return abs((1/2)*A) def new_pgon(self, current_line): """Takes as input a line and creates a new polygon from the intersects. Outputs a list.""" Q=[] if current_line.b==0: for i in range(len(self.vertex_list)): if i < len(self.vertex_list)-1: #does_intersect returns intersection if it exists r = self.does_intersect(self.vertex_list[i],self.vertex_list[i+1],current_line) if r != None: # meaning there exists intersection Q.append(r) g=self.vertex_list[i+1] if g[0] <= current_line.c: Q.append(g) if i == len(self.vertex_list)-1: r= self.does_intersect(self.vertex_list[i],self.vertex_list[0],current_line) if r!= None: Q.append(r) g=self.vertex_list[0] if g[1] < current_line.c: Q.append(g) N = pgon(Q) return N else: for i in range(len(self.vertex_list)): if i < len(self.vertex_list)-1: #does_intersect returns intersection if it exists r = self.does_intersect(self.vertex_list[i+1],self.vertex_list[i],current_line) if r != None: # meaning there exists intersection Q.append(r) g=self.vertex_list[i+1] if g[1] < -(current_line.a/current_line.b)*g[0]+(current_line.c/current_line.b): Q.append(g) if i == len(self.vertex_list)-1: r = self.does_intersect(self.vertex_list[0],self.vertex_list[i],current_line) if r!= None: Q.append(r) g=self.vertex_list[0] if g[1] < -(current_line.a/current_line.b)*g[0]+(current_line.c/current_line.b): Q.append(g) N = pgon(Q) return N def cent_x(self): """Takes a list of Pgon vertices and finds the x-centroid""" S=[] for i in range(0,len(self.vertex_list)): for z in self.vertex_list[i]: S.append(z) B=self.area if B==0: return A=0 for i in range(0,len(S),2): if i < len(S)-2: A= A + (S[i]+S[i+2])*(S[i]*S[i+3] - S[i+2]*S[i+1]) else: A= A + (S[i]+S[0])*(S[i]*S[1] - S[0]*S[i+1]) return (1/(6*B))*A def cent_y(self): """Takes a list of polygon vertices and finds the y-centroid""" S=[] for i in range(0,len(self.vertex_list)): for z in self.vertex_list[i]: S.append(z) B=self.area if B==0: return A=0 for i in range(0,len(S),2): if i < len(S)-2: A= A + (S[i+1]+S[i+3])*(S[i]*S[i+3] - S[i+2]*S[i+1]) else: A= A + (S[i+1]+S[1])*(S[i]*S[1] - S[0]*S[i+1]) return (1/(6*B))*A def does_intersect(self, pt1,pt2,current_line): """Determines if two lines intersect. Returns point of intersection.""" line1=line(pt1,pt2).give() ###.give() returns a list! [a,b,c] line2=current_line.give() m=matrix([line1,line2]) x0=m.solve_right(m)[0,2] y0=m.solve_right(m)[1,2] if pt1[0]!=pt2[0]: if pt1[0]>=pt2[0]: if pt2[0] <= x0 and pt1[0] >= x0: return (x0,y0) if pt1[0]<= pt2[0]: if pt1[0] <= x0 and pt2[0] >= x0: return (x0,y0) else: if pt1[1] <= pt2[1]: if pt2[1] >= y0 and pt1[1] <= y0: return (x0,y0) if pt1[1] >= pt2[1]: if pt1[1] >= y0 and pt2[1] <= y0: return (x0,y0) class line: def __init__(self, p1 = (0,0), p2 = (0,0) ): if p2[0]-p1[0]==0: self.a=1 self.b=0 self.c=p1[0] else: d = ( (p2[1] - p1[1]) / (p2[0] - p1[0]) ) self.a = -d self.b = 1 self.c = p1[1] - (d*p1[0]) def give(self): ls = [self.a, self.b, self.c] return ls class hamline: def __init__(self, xmin, xmax, ymin, ymax, resolution, polygons, HALF_AREA): self.steps = int(1/resolution) #Pass in resolution as a percent (meaning less than 1)! self.resolution = resolution self.xmin = xmin self.xmax = xmax self.ymin = ymin self.ymax = ymax self.polygons = polygons #pass this in a list of two pgons self.HALF_AREA = HALF_AREA self.found = False self.areaLine = self.find_line() self.topSearch = False def __repr__(self): print "%sy + %sx = %s"%(self.a, self.b, self.c) return "[%s, %s, %s]"%(self.a, self.b, self.c) ###This returns a line: def find_line(self): steps = self.steps xmax = self.xmax xmin = self.xmin ymax = self.ymax ymin = self.ymin ###These are all just points stored as tuples corner1 = (xmax, ymax) corner2 = (xmax, ymin) corner3 = (xmin, ymin) corner4 = (xmin, ymax) divis = [0..steps] comps = [] #as in components yrange = ymax-ymin xrange = xmax-xmin ###The thing is that all of our lines that we are going to use will be generated ###Generate points along the line ### leftline_points etc are lists of 2-tuples leftline_points = [] rightline_points = [] topline_points = [] bottomline_points = [] ####Note: As implemented, we divide the region into step by step components. #### In a case where xrange >> yrange or vice versa, this could be a problem #### I suggest doing component wise entry of steps in the input box for i in divis: leftline_points.append( (xmin, (i*(yrange/steps) + ymin)) ) rightline_points.append( (xmax, (i*(yrange/steps) + ymin)) ) topline_points.append( ( (i*(xrange/steps) + xmin), ymax) ) bottomline_points.append( ( (i*(xrange/steps) + xmin), ymin) ) #### Post-cond: we have all our points from which to make our lovely lines! #### Check the lines from leftline to rightline for i in leftline_points: for j in rightline_points: connect = line(i,j) #We have the set of two pgons. Those are invariant. New pgons get generated. Find areas of those: # Each pgon has a vertex_list that can create a new pgon # newpgon takes a vertex_list and a line Poly1 = self.polygons[0].vertex_list Poly2 = self.polygons[1].vertex_list P1 = pgon(Poly1) P2 = pgon(Poly2) P1_intersect = P1.new_pgon(connect) ## these are new polygons, created by line P2_intersect = P2.new_pgon(connect) ALP1= P1_intersect.area AP1 = P1.area ALP2 = P2_intersect.area AP2 = P2.area #This is assuming that the new pgons both appear below the line if not( (P1_intersect.area==0) or (P2_intersect.area ==0) ) and not( (P1_intersect.area == P1.area) or (P2_intersect.area == P2.area) ): if ( (ALP1/AP1 < self.resolution + 1/2) and (ALP1/AP1 > 1/2 - self.resolution) ): if ( (ALP2/AP2 < self.resolution + 1/2 ) and (ALP2/AP2 > - self.resolution + 1/2) ): self.found = True return connect if (not self.found): self.topSearch = True for i in topline_points: for j in bottomline_points: connect = line(i,j) #We have the set of two pgons. Those are invariant. New pgons get generated. #Find areas of those: # Each pgon has a vertex_list that can create a new pgon # newpgon takes a vertex_list and a line Poly1 = self.polygons[0].vertex_list Poly2 = self.polygons[1].vertex_list P1 = pgon(Poly1) P2 = pgon(Poly2) P1_intersect = P1.new_pgon(connect) P2_intersect = P2.new_pgon(connect) ALP1= P1_intersect.area AP1 = P1.area ALP2 = P2_intersect.area AP2 = P2.area #This is assuming that the new pgons both appear below the line if not( (P1_intersect.area==0) or (P2_intersect.area ==0) ) and not( (P1_intersect.area == P1.area) or (P2_intersect.area == P2.area) ): if ( (ALP1/AP1 < self.resolution + 1/2) and (ALP1/AP1 > 1/2 - self.resolution) ): if ( (ALP2/AP2 < self.resolution + 1/2 ) and (ALP2/AP2 > - self.resolution + 1/2) ): self.found = True return connect if (not self.found): Poly1 = self.polygons[0].vertex_list Poly2 = self.polygons[1].vertex_list P1 = pgon(Poly1) P2 = pgon(Poly2) C1= P1.centroid C2 = P2.centroid return line(C1, C2) @interact def f(vertices=input_box("(1,-2) (2,2) (4,-1) (6,5) (-1,2)", "Vertices P1 (separated by spaces)", type=str), clr=color_selector(Color('purple'), widget='jpicker'), vertices2=input_box("(-1,-1) (-2,-2) (-5,-1/2) (-5,-8) (-4,-3) (-3,-7) (-5/2,-2) (-2,-6) (0,-1)", "Vertices P2 (separated by spaces)", type=str), clr2=color_selector(Color('red'), widget='jpicker'), error=slider(0,.5,default=.2)): v = sage_eval('['+vertices.replace(')','),')+']') v2 = sage_eval('['+vertices2.replace(')','),')+']') """Debugging""" #P1Plot = polygon(v, color=clr) #P2Plot = polygon(v2, color=clr2) #show(P1Plot + P2Plot) #### Storing the information from the vertices inputed into two shapes. shape1 = pgon(v) shape2 = pgon(v2) ####Areas of initial polygons. A1= shape1.area Half_A1=A1/2 A2= shape2.area Half_A2= A2/2 #### Centroids of polygons and connecting line. C1= shape1.centroid C2= shape2.centroid CLine = line(C1, C2) L= CLine m1= max_list(v,v2) n1=min_list(v,v2) ####Form a new polygon from intersecting original with centroid line. CP1=shape1.new_pgon(L) CP2=shape2.new_pgon(L) ####Areas of bottom polygons obtained from the centroid line. CA1=float(CP1.PArea()) CA2=float(CP2.PArea()) bottom = CA1 + CA2 top = (A1-CA1) + (A2-CA2) HALF_AREA = (A1 +A2)/2.0 resolution = error ####Calculate the x and y max of the polygon vertice sets. To be used for constructing the plot. ymax = find_ymax(v, v2) ymin = find_ymin(v, v2) xmax = find_xmax(v, v2) xmin = find_xmin(v, v2) ####Create hamline, then given hamline defines the new polygons below the hamline. poly_list = [shape1, shape2] ham = hamline(xmin, xmax, ymin, ymax, resolution, poly_list, HALF_AREA) areaCutLine = ham.areaLine if (areaCutLine) != None: L2 = areaCutLine.give() else: P1Plot = polygon(v, color=clr) P2Plot = polygon(v2, color=clr2) CLinePlot = plot(-(L.a/L.b)*x+L.c/L.b,n1-1,m1+1, thickness=.3, color='black') show(P1Plot + P2Plot + CLinePlot) print 'Self-intersecting polygons not allowed.' return #import sys #sys.exit ("Your Polygon Has a Self Intersection.") NP1=shape1.new_pgon(ham.areaLine) NP2=shape2.new_pgon(ham.areaLine) ####Areas of the newly formed bottom polygons. BA1=float(NP1.PArea()) BA2=float(NP2.PArea()) bottom = BA1 + BA2 top = (A1-BA1) + (A2-BA2) HALF_AREA = (bottom + top)/2.0 ####Polygons and lines to be plotted. P1Plot = polygon(v, color=clr) P2Plot = polygon(v2, color=clr2) CLinePlot = plot(-(L.a/L.b)*x+L.c/L.b,n1,m1, thickness=.3, color='black') if -(L2[0]/L2[1])*xmax+L2[2]/L2[1] > ymax: HAMLINEPlot= plot(-(L2[0]/L2[1])*x+L2[2]/L2[1],( ymin - (L2[2]/L2[1]) ) / (-L2[0]/L2[1] ),( ymax - (L2[2]/L2[1]) ) / (-L2[0]/L2[1] ), thickness=2, color='black') else: HAMLINEPlot= plot(-(L2[0]/L2[1])*x+L2[2]/L2[1],n1,m1, thickness=2, color='black') if NP1.vertex_list: P1BPLOT = polygon(NP1.vertex_list, color= 'green') else: print 'No Bottom Found for P1, check for self-intersection or change error.' show(P1Plot+ P2Plot) return if NP2.vertex_list: P2BPLOT = polygon(NP2.vertex_list, color= 'yellow') else: print 'No Bottom Found for P2, check for self-intersection or change error.' show(P1Plot+ P2Plot) return #### Displayed Graph and Quality of Approximation Printouts # show( 'The line: y=%s results in P1 error= %s and P2 error = %s'%((-(L2[0]/L2[1])*x+L2[2]/L2[1]),(round(.5 -(BA1)/A1, 3)),(round(.5-(BA2)/A2, 3))) ) show(HAMLINEPlot + P1Plot+ P2Plot +P1BPLOT + P2BPLOT) 
       

Click to the left again to hide and once more to show the dynamic interactive window

# una cadena de Markov html('<h2>El casino que (casi) no hace trampa</h2>') mat = matrix([[0.8,0.2],[0.2,0.8]]) # cargado / legal m = hmm.DiscreteHiddenMarkovModel(mat, [[1/7,1/7,1/7,1/7,1/7,2/7],[1/6,1/6,1/6,1/6,1/6,1/6]], [.2,.8],emission_symbols=[1,2,3,4,5,6]) m = hmm.DiscreteHiddenMarkovModel([[0.75,0.25],[0.2,0.8]], [[1/10,1/10,1/10,1/10,1/10,1/2],[1/6,1/6,1/6,1/6,1/6,1/6]], [.2,.8],emission_symbols=[1,2,3,4,5,6]) #show(m.graph().plot(), figsize=4) sample = m.generate_sequence(10) print "Resultados (visibles): ", sample[0] print "Tipo de dados (oculto): ", sample[1] 
       

A continuación veremos:

  1. Los resultados de una serie de tiradas según las reglas anteriores.
  2. Qué dado se usó en cada tirada (uno equilibrado o uno cargado).
  3. Qué predice el algoritmo de Viterbi a partir de los resultados visibles.
@interact def dishonest_casino(auto_update=False): test = list(m.generate_sequence(75)) vit_test = list(m.viterbi(test[0])[0]) html('<h3>El casino que (casi) no hace trampa</h3>') html('Tiradas: '+str(test[0]).replace(',','').replace(' ','')[1:-1]) vit_str = str(vit_test).replace(',','').replace(' ','') vit_str = vit_str.replace('1','F').replace('0','<font color="#FF0000">L</font>')[1:-1] actual_str = str(list(test[1])).replace(',','').replace(' ','') actual_str = actual_str.replace('1','F').replace('0','<font color="#FF0000">L</font>')[1:-1] html(' Real: '+ actual_str) html('Viterbi: '+vit_str) 
       

Click to the left again to hide and once more to show the dynamic interactive window