{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lista 2 Física computacional "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "import matplotlib.pyplot as plt\n",
    "# suporte para graficos 3D\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "plt.rc('text', usetex=True)\n",
    "plt.rcParams.update({'font.size': 18})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mapa logístico\n",
    "\n",
    "\n",
    "Objetivos de hoje:\n",
    "\n",
    "- *familiarização com sistemas dinâmicos*\n",
    "    - iteração    \n",
    "    - atratores\n",
    "    - períodos e ciclos\n",
    "    - caos\n",
    "- *mapa logístico*\n",
    "\n",
    "$$ x_{k+1} = \\rho \\, x_k (1 - x_k)  $$\n",
    "\n",
    "para $k=0,1,\\ldots$, com $x \\in [0,1]$ e $\\rho \\in [0,4] $.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercício\n",
    "\n",
    "1. Descreva brevemente o mapa logístico. \n",
    "2. Quais os tipos de problema podem ser modelados pelo mapa logístico?\n",
    "3. Exemplifique."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "escreva sua resposta aqui\n",
    "\n",
    "----\n",
    "\n",
    "- 1.\n",
    "\n",
    "----\n",
    "\n",
    "- 2.\n",
    "\n",
    "----\n",
    "\n",
    "- 3.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercício\n",
    "\n",
    "Comente os códigos computacionais abaixo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mapa_logistico(x,params):\n",
    "    \"\"\"\n",
    "    retorna uma única interação do mapa logistico.\n",
    "    \n",
    "    x     : real entre 0 e 1\n",
    "    params: **lista** contendo um único número real entre 0 e 4\n",
    "    \"\"\"\n",
    "    return params[0]*x*(1e0-x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def iteracao(f,N=100, x0=0.123456 ,params=[3.0]):\n",
    "    \"\"\"\n",
    "    retorna a N-ésima iteração do mapa f, partindo da semente x0, com a lista de params\n",
    "    \n",
    "    f     : nome da função que descreve o mapa\n",
    "    N     : número de iterações do sistema dinâmico\n",
    "    x0    : semente\n",
    "    params: lista com parâmetros do mapa f\n",
    "    \n",
    "    y     : saída. 1D-array de tamanho N\n",
    "    \"\"\"\n",
    "    y = np.zeros(N)\n",
    "    y[0] = x0\n",
    "    for k in range(N-1):\n",
    "        y[k+1] = f(y[k],params)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercício \n",
    "\n",
    "A partir da função `iteracao`, defina seu método de iteração que guarde apenas as $m$ últimas iterações."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def iteracao_otimizada(f,N=100,x0=0.123456, params=[1.0],m=128):\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercício \n",
    "\n",
    "1. Comente o código computacional abaixo\n",
    "2. Qual o objetivo dos gráficos abaixo?\n",
    "3. Existe um ponto fixo?\n",
    "4. Varie a condição inicial $x_0$. O ponto fixo muda?\n",
    "5. Repita os passos 3 e 4 para $\\rho=3.25$ e $3.75$ com $N=2000$ iterações."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "total = 100\n",
    "t = np.arange(total)\n",
    "x = iteracao(mapa_logistico,x0=0.23,params=[0.8],N=total);\n",
    "\n",
    "plt.figure(figsize=(18,6))\n",
    "plt.subplot(1,2,1)\n",
    "plt.xlabel('$t$')\n",
    "plt.ylabel('$x(t)$')\n",
    "plt.scatter(t,x,c=t); \n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "h,b = np.histogram(x,range=[0,1],bins=50);\n",
    "plt.bar(b[:-1],h, width = 0.01) ; \n",
    "plt.xlabel('$x$')\n",
    "plt.ylabel('frequencia')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rho 3.25"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rho 3.75"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercícios\n",
    "\n",
    "1. Comente o código abaixo.\n",
    "2. Construa o diagrama de pontos fixos para o mapa logístico."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def diagrama(mapa,amostras=128,delta=0.01,params=[0.0,4.0],x0=0.49,N=2000):\n",
    "    \"\"\"\n",
    "    coloque a descrição geral das funções aqui. Essa região é o __doc__ da função, \n",
    "    que contém a documentação, incluindo objetivo da função, parâmetros e retornos.\n",
    "    \"\"\"\n",
    "    amostras = 128\n",
    "    plt.figure(figsize=(12,8))\n",
    "    plt.grid(color='black', linestyle='-', linewidth=1,alpha=0.25)\n",
    "    for rho in  np.arange(params[0],params[1],delta):\n",
    "        x = iteracao(mapa,x0=x0,params=[rho],N=N)[-amostras:];\n",
    "        plt.scatter(rho*np.ones(amostras), x,c='purple',s=15,marker='.',alpha=0.75);\n",
    "    plt.ylabel(u'$x^*$')    \n",
    "    plt.xlabel(u'$\\\\rho$')\n",
    "    return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# grafico do ponto fixo contra diferentes valores do parametro\n",
    "diagrama(mapa_logistico,N=2000,amostras=128);\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercício\n",
    "\n",
    "1. Leia e comente o código abaixo\n",
    "2. Interprete a saída.\n",
    "3. Execute a função `gera_relatorio` para $\\rho = 1.5, 3.1, 3.5$ e $3.9$. Faça uma breve descrição de cada caso.\n",
    "4. Modifique a função `gera_relatorio` de maneira que ela também retorne a entropia de Boltzmann-Shannon, $$S = -\\sum_{k} p_k \\ln p_k, $$ em que $k$ especifica o intervalo de classe (bin) e $p_k$ corresponde à probabilidade de observar a classe $k$. Quais valores de $\\rho$ apresentam maior entropia?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gera_relatorio(mapa,params,x0=0.99,N=1000):\n",
    "    x = iteracao(mapa,x0=x0,params=params,N=N);\n",
    "    plt.rcParams.update({'font.size': 18})    \n",
    "    plt.subplots(nrows=2,ncols=2,figsize=(18,12))    \n",
    "    plt.tight_layout() \n",
    "    \n",
    "    plt.subplot(2,2,1)    \n",
    "    plt.scatter(np.arange(len(x)),x,c=np.arange(len(x)))\n",
    "    plt.title('$x_t$')\n",
    "    plt.xlabel('$t$')\n",
    "    plt.ylabel('$x_{t}$')\n",
    "    plt.legend(['$\\\\rho = %s $' % (str(params))])\n",
    "            \n",
    "    plt.subplot(2,2,2)    \n",
    "    plt.scatter(x[:-1], x[1:],marker='o',s=20*np.log(np.arange(1,N))+1,alpha=0.5);\n",
    "    plt.title('$x_t \\\\times x_{t+1}  $')\n",
    "    plt.legend(['$\\\\rho = %s $' % (str(params))])\n",
    "    plt.xlabel('$x_{t}$')\n",
    "    plt.ylabel('$x_{t+1}$')\n",
    "\n",
    "    ax = plt.subplot(2,2,3,projection='3d')\n",
    "    ax.scatter(x[:-2],x[1:-1],x[2:],s=100,c=np.arange(0,len(x[:-2])))\n",
    "    angle=20\n",
    "    ax.view_init(30, angle)\n",
    "    ax.scatter(x[:-2],x[1:-1],c='b',alpha=0.25)\n",
    "    \n",
    "    plt.subplot(2,2,4)    \n",
    "    h,b = np.histogram(x,range=[0,1],bins=50);\n",
    "    plt.bar(b[1:]*.5+b[:-1]*.5,h, width = 0.01)\n",
    "    plt.title('histograma')\n",
    "    plt.legend(['$\\\\rho = %s $' % (str(params))])\n",
    "    plt.ylabel('frequencia')\n",
    "    plt.xlabel('$x$')\n",
    "    \n",
    "    plt.tight_layout()     \n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rho = 1.5\n",
    "x = gera_relatorio(mapa_logistico,[1.5],N=1000)    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rho = 3.1\n",
    "x = gera_relatorio(mapa_logistico,[3.1],N=1000)    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rho = 3.5\n",
    "x = gera_relatorio(mapa_logistico,[3.5],N=1000)    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rho = 3.9\n",
    "x = gera_relatorio(mapa_logistico,[3.9],N=5000)    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercício\n",
    "\n",
    "Considere o mapa $$x_{k+1} =  x_k \\,\\textrm{e}^{-\\rho(1- x_k)},$$\n",
    "em que $\\rho\\geqslant 0$ e $x \\geqslant 0$.\n",
    "\n",
    "1. Qual a interpretação desse mapa? \n",
    "2. Escreva o código computacional do mapa acima.\n",
    "3. Use a função `gera_relatorio` para $x_0=0.01$, $\\rho = 0.01$, com $N=1000$. Repita o procedimento com $x_0=0.99$. O ponto fixo é atrativo, repulsivo ou marginal?\n",
    "4. Construa o diagrama do mapa ( pontos fixos $x^*$ por $\\rho$ ). Utilize o incremento $\\Delta \\rho = 0.01$ para $\\rho \\in [0, 5]$ e use $m=128$ amostras.\n",
    "5. Calcule ou estime os pontos fixos para $\\rho =5/2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. coloque sua interpretação aqui"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ex 2\n",
    "def mapa_ricker(x,params=[1.0]):\n",
    "    pass    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ex 3\n",
    "x = gera_relatorio(mapa_ricker,[0.01],N=1000,x0=0.01)    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ex 4\n",
    "diagrama(mapa_ricker,params=[0,4],x0=3.01);\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ex 5"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}