regressões com unidades

LabIFSC2/_arrays.py

@@ -31,8 +31,10 @@ def nominais(arrayMedidas : np.ndarray,unidade:str) -> np.ndarray:
     '''
     if not (isinstance(arrayMedidas[0],Medida)):
         raise TypeError('Os valores do array não são Medidas')
-
-    return np.array([medida._nominal.to(unidade).magnitude for medida in arrayMedidas],dtype=float)
+    if unidade=='si':
+        return np.array([medida._nominal.to_base_units().magnitude for medida in arrayMedidas],dtype=float)
+    else:
+        return np.array([medida._nominal.to(unidade).magnitude for medida in arrayMedidas],dtype=float)
 
 @obrigar_tipos
 def incertezas(arrayMedidas : np.ndarray,unidade:str) -> np.ndarray:
@@ -54,7 +56,10 @@ def incertezas(arrayMedidas : np.ndarray,unidade:str) -> np.ndarray:
 '''
     if not (isinstance(arrayMedidas[0],Medida)):
         raise TypeError('Os valores do array não são Medidas')
-    return np.array([medida.incerteza(unidade) for medida in arrayMedidas],dtype=float)
+    if unidade=='si':
+        return np.array([medida._incerteza.to_base_units().magnitude for medida in arrayMedidas],dtype=float)
+    else:
+        return np.array([medida._incerteza.to(unidade).magnitude for medida in arrayMedidas],dtype=float)
 
 @obrigar_tipos
 def curva_min(arrayMedidas : np.ndarray,unidade:str,sigma: float | int =2) -> np.ndarray:

LabIFSC2/_regressões.py

@@ -1,9 +1,10 @@
 import string
-from collections.abc import Iterator
+from collections.abc import Callable, Iterator
 from numbers import Real
 from typing import Any
 
 import numpy as np
+import pint
 from numpy.polynomial import Polynomial
 
 from ._arrays import arrayM, nominais
@@ -12,6 +13,28 @@
 from ._tipagem_forte import obrigar_tipos
 
 
+def _aplicar_funcao_sem_passar_pelo_sistema_de_unidades(
+        array_medidas:np.ndarray,lab_func:Callable)->np.ndarray:
+    '''
+    #precisamos aplicar log/exp sem passar pelo sistema de unidades
+    #um método meio quicky and dirty mas necessário para não dar erro de dimensão
+    '''
+    medidas_novas=[]
+    for medida in array_medidas:
+        unidade=str(medida._nominal.units)
+        medida_intermediaria=Medida(medida._nominal.magnitude,medida._incerteza._magnitude,"")
+        medida_intermediaria=lab_func(medida_intermediaria)
+        nominal=medida_intermediaria._nominal._magnitude
+        incerteza=medida_intermediaria._incerteza._magnitude
+        medidas_novas.append(Medida(nominal,incerteza,unidade))
+    return np.array(medidas_novas)
+def _forcar_troca_de_unidade(medida:Medida|np.ndarray,unidade:str)->Medida|np.ndarray:
+    if isinstance(medida,np.ndarray):
+        return np.array([Medida(med._nominal.magnitude,med._incerteza.magnitude,unidade) for med in medida])
+    else:
+        return Medida(medida._nominal.magnitude,medida._incerteza.magnitude,unidade)
+
+
 class MPolinomio:
     @obrigar_tipos
     def __init__(self,coeficientes:np.ndarray):
@@ -51,7 +74,7 @@ def __init__(self,a:Medida,k:Medida,base:Real):
 
     @obrigar_tipos
     def __call__(self,x:Medida | np.ndarray) -> Medida | np.ndarray:
-        resultado:  Medida | np.ndarray=power(float(self.base),x*self.k)
+        resultado:  Medida | np.ndarray=self.a*power(float(self.base),x*self.k)
         return resultado
 
     def __repr__(self)->str:
@@ -63,25 +86,39 @@ class MLeiDePotencia:
     '''Classe para modelar uma função de lei de potência
     y = a * x^n
     '''
-    __slots__ = ['a', 'n','_valores']
 
     @obrigar_tipos
-    def __init__(self, a: Medida, n: Medida):
+    def __init__(self, a: Medida, n: Medida,y_unidade:pint.Quantity):
         self.a = a
         self.n = n
         self._valores = (a, n)
+        self._y_unidade=y_unidade
 
     @obrigar_tipos
     def __call__(self, x:Medida | np.ndarray) -> Medida | np.ndarray:
-        resultado : Medida | np.ndarray=self.a *power(x, self.n)
+        #Aqui teremos um pequeno passo quick and dirty de novo
+        #se tentarmos fazer a exponencial de um PlainQuatity temos um erro:
+        #PlainQuantity array exponents are only allowed if the plain is dimensionless
+        #então precisamos fazer umas conversões manuais aqui
+        if isinstance(x,np.ndarray): unidade=(x[0]._nominal**self.n._nominal)
+        else: unidade=(x._nominal**self.n._nominal)
+        x=_forcar_troca_de_unidade(x,'dimensionless')
+        parte_exponencial=x**self.n
+        resultado:Medida | np.ndarray=self.a*_forcar_troca_de_unidade(parte_exponencial,str(unidade.units))
+
+        resultado_pint=resultado._nominal if isinstance(resultado,Medida) else resultado[0]._nominal
+        if not resultado_pint.is_compatible_with(self._y_unidade):
+            raise ValueError("Unidades incompatíveis")
         return resultado
 
     def __repr__(self)->str:
         return f'MLeiDePotencia(a={self.a}, b={self.n})'
 
     def __iter__(self)->Iterator[object]:
         return iter(self._valores)
-
+
+
+
 @obrigar_tipos
 def regressao_polinomial(x_medidas:np.ndarray,y_medidas:np.ndarray,grau:int) -> MPolinomio:
     if len(x_medidas)!=len(y_medidas):
@@ -90,11 +127,15 @@ def regressao_polinomial(x_medidas:np.ndarray,y_medidas:np.ndarray,grau:int) ->
         raise ValueError("Não há dados suficientes para um polinômio de grau {grau} (overfitting)")
     if not (isinstance(x_medidas[0],Medida) and isinstance(y_medidas[0],Medida)):
         raise TypeError('x_medidas e y_medidas precisam ser np.ndarray de medidas')
-    x_medidas=np.array([x._nominal.magnitude for x in x_medidas])
-    y_medidas=np.array([y._nominal.magnitude for y in y_medidas])
-    p, cov = np.polyfit(x_medidas.astype(float), y_medidas.astype(float), grau, cov=True)
-    medidas_coeficientes = np.array([Medida(valor, np.sqrt(cov[i, i]),'') for i, valor in enumerate(p)],dtype=Medida)
-    return MPolinomio(medidas_coeficientes)
+    x_float=np.array([x._nominal.magnitude for x in x_medidas],dtype=float)
+    y_float=np.array([y._nominal.magnitude for y in y_medidas],dtype=float)
+    p, cov = np.polyfit(x_float, y_float, grau, cov=True)
+    erros=np.sqrt(np.diag(cov))
+    medidas_coeficientes=[]
+    for index in range(len(p)):
+        unidade= str((y_medidas[0]._nominal/x_medidas[0]._nominal**(grau-index)).units)
+        medidas_coeficientes.append(Medida(p[index],erros[index],unidade))
+    return MPolinomio(np.array(medidas_coeficientes))
 
 @obrigar_tipos
 def regressao_linear(x_medidas:np.ndarray,
@@ -111,17 +152,26 @@ def regressao_exponencial(x_medidas:np.ndarray,y_medidas:np.ndarray,
 
     if base<1: raise ValueError('Base precisa ser maior que 1')
 
-    polinomio=regressao_linear(x_medidas,log(y_medidas)/log(base))
+    pegar_log=lambda x: log(x)/log(base)
+    log_y_medidas=_aplicar_funcao_sem_passar_pelo_sistema_de_unidades(y_medidas,pegar_log) 
+    polinomio=regressao_linear(x_medidas,log_y_medidas)
     k=polinomio.a
-    a=exp(polinomio.b)
+    a=_aplicar_funcao_sem_passar_pelo_sistema_de_unidades(np.array([polinomio.b]),exp)[0]
+
+    k=_forcar_troca_de_unidade(k,str((1/x_medidas[0]._nominal).units))
+    a=_forcar_troca_de_unidade(a,str(y_medidas[0]._nominal.units))
     return MExponencial(a,k,base)
 
 @obrigar_tipos
 def regressao_potencia(x_medidas:np.ndarray, y_medidas:np.ndarray) -> MLeiDePotencia:
 
     if not bool(np.all([y._nominal.magnitude>0 for y in y_medidas]) and np.all([x._nominal.magnitude>0 for x in x_medidas])):
             raise ValueError('Todos x e y precisam ser positivos para uma modelagem exponencial')
-    polinomio=regressao_linear(log(x_medidas),log(y_medidas))
-    a=exp(polinomio.b)
+    log_y_medidas=_aplicar_funcao_sem_passar_pelo_sistema_de_unidades(y_medidas,log)
+    log_x_medidas=_aplicar_funcao_sem_passar_pelo_sistema_de_unidades(x_medidas,log)
+    polinomio=regressao_linear(log_x_medidas,log_y_medidas)
+    a=_aplicar_funcao_sem_passar_pelo_sistema_de_unidades(np.array([polinomio.b]),exp)[0]
     n=polinomio.a
-    return MLeiDePotencia(a,n)
+    a=_forcar_troca_de_unidade(a,str((y_medidas[0]._nominal/x_medidas[0]._nominal**n._nominal.magnitude).units))
+    n=_forcar_troca_de_unidade(n,"dimensionless")
+    return MLeiDePotencia(a,n,y_medidas[0]._nominal)

docs/Arrays/get_incertezas_nominais.md

@@ -1,4 +1,4 @@
 
-:::LabIFSC2._operacoes_em_arrays.incertezas
-:::LabIFSC2._operacoes_em_arrays.nominais
+:::LabIFSC2._arrays.incertezas
+:::LabIFSC2._arrays.nominais
 

docs/Arrays/linspace.md

@@ -1,2 +1,2 @@
-:::LabIFSC2._operacoes_em_arrays.linspace
+:::LabIFSC2._arrays.linspace
 

docs/graficos.md

@@ -2,7 +2,27 @@ O LabIFSC2 precisa implementar algumas funções para converter um objeto `Medid
 float, para que se possa criar gráficos em bibliotecas como o [Matplotlib](https://atplotlib.org/)/[Seaborn](https://seaborn.pydata.org/). Nessa secção continuaremos o exemplo do campo magnético em função da distância da secção [Arrays](arrays.md)
 
 ## Nominais
+Para obter os valores nominais de um array numpy de medidas basta usar a função `nominais(array_medida,unidade)`
+```py
+--8<-- "tests/test_doc_nominal.py:1:5"
+```
 
 
+## Incertezas
+De maneira análogo podemos também pegar as incertezas com `incertezas(array_medida,unidade)`
+```py
+--8<-- "tests/test_doc_incerteza.py:1:5"
+```
 
-## Incertezas
+## Dispersão com barras de erro 
+Eis um exemplo simples de como fazer um gráfico de dispersão com erros tanto em x quanto em y,
+basta usar a função do matplotlib [`plt.errorbar`](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.errorbar.html), `nominais` e `incertezas`.
+
+```py
+--8<-- "tests/test_doc_grafico_scatter.py:1:16"
+```
+
+Repare que as unidade são variáveis no código que podem ser modificadas rapidamente, o resultado é esse:
+<img src="exemplo.jpg" width=400>
+
+##  

docs/introducao.md

@@ -29,6 +29,14 @@ Com o LabIFSC2, podemos fazer os cálculos assim:
 ```
 Calculamos então que o IMC tem valor esperado de \(\mu=24,5 \, kg/m²\) e desvio padrão de \(\sigma=0.3 \, kg/m²\).
 
+Você pode acessar o valor nominal e a incerteza fazendo `medida.nominal(unidade)` ou `medida.incerteza(unidade)`,
+repare que é necessário especificar uma unidade, visto que o valor nominal/incerteza são totalmente dependentes
+dela. É importante notar que somente em casos bem específicos você vai acessar diretamente esses valores, se você usa muito essas funções talvez não esteja usando corretamente a biblioteca
+
+```py title="Cálculo de IMC"
+--8<-- "tests/test_doc_imc.py:7:8"
+```
+
 ## Comparando Medidas
 Se uma segunda pessoa afirmar que seu IMC é de (25 ± 0.1), podemos dizer que seus IMCs são equivalentes? Mesmo que \(25 \ne 24.5\), pela incerteza nas medidas podemos dizer que essa diferença está na margem de erro do experimento.
 

pyproject.toml

@@ -34,6 +34,10 @@ taskipy = "^1.14.1"
 isort = "^5.13.2"
 mypy = "1.14.1"
 
+
+[tool.poetry.group.testes.dependencies]
+matplotlib = "^3.10.0"
+
 [tool.mypy]
 strict = true
 disallow_any_generics = false

tests/test_arrayM.py

@@ -26,6 +26,7 @@ def test_incerteza_array():
     assert incerteza[3]==0
     with pytest.raises(TypeError):
         incertezas(np.arange(10),'')
+
 
 def test_curvamin():
     t=np.array([Medida(5,0.1,''),Medida(9,2,''),Medida(11,0.5,'')])

tests/test_doc_grafico_fitting.py

@@ -0,0 +1,25 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+from LabIFSC2 import *
+
+campo_magnético=arrayM([210,90,70,54,39,32,33,27,22,20],1,'muT')
+distancias=linspace(1,10,10,0.01,'cm') 
+
+unidade_x='cm'
+unidade_y='muT'
+
+regressao=regressao_potencia(distancias,campo_magnético)
+print(regressao)
+plt.style.use('ggplot')
+plt.errorbar(nominais(distancias,unidade_x),nominais(campo_magnético,unidade_y),
+             xerr=incertezas(distancias,unidade_x),yerr=incertezas(campo_magnético,unidade_y),
+             fmt='o',label='Dados experimentais',color='red')
+
+x=linspace(1,10,100,0,unidade_x)
+plt.plot(nominais(x,unidade_x),nominais(regressao(x),unidade_y),color='blue',
+         label="Curva teórica")
+plt.fill_between(nominais(x,unidade_x),curva_min(regressao(x),unidade_y),
+                 curva_max(regressao(x),unidade_y),color='blue',alpha=0.3)
+plt.legend()
+plt.savefig('teste.jpg') 

tests/test_doc_grafico_scatter.py

@@ -0,0 +1,17 @@
+import matplotlib.pyplot as plt
+
+from LabIFSC2 import *
+
+campo_magnético=arrayM([250,150,110,90,70,60,55,40,25,20],1,'muT')
+distancias=linspace(1,10,10,0.5,'cm') 
+
+unidade_x='cm'
+unidade_y='muT'
+
+plt.style.use('ggplot')
+plt.errorbar(nominais(distancias,unidade_x),nominais(campo_magnético,unidade_y),
+             xerr=incertezas(distancias,unidade_x),yerr=incertezas(campo_magnético,unidade_y),fmt='o')
+
+plt.xlabel(f"Distancia ({unidade_x})")
+plt.ylabel(f"Campo magnético ({unidade_y})")
+#plt.savefig("exemplo.jpg")

tests/test_doc_incerteza.py

@@ -0,0 +1,7 @@
+from LabIFSC2 import *
+
+campo_magnético=arrayM([250,150,110,90,70,60,55,40,25,20],1,'muT')
+print(incertezas(campo_magnético,'muT'))
+#[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
+
+assert str(incertezas(campo_magnético,'muT'))=="[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]"

tests/test_doc_nominal.py

@@ -0,0 +1,7 @@
+from LabIFSC2 import *
+
+campo_magnético=arrayM([250,150,110,90,70,60,55,40,25,20],1,'muT')
+print(nominais(campo_magnético,'muT'))
+#[250. 150. 110.  90.  70.  60.  55.  40.  25.  20.]
+
+assert str(nominais(campo_magnético,'muT'))=="[250. 150. 110.  90.  70.  60.  55.  40.  25.  20.]"

tests/test_regressao_unidades.py

@@ -0,0 +1,39 @@
+import pint
+import pytest
+
+from LabIFSC2 import *
+
+campo_magnético=arrayM([210,90,70,54,39,32,33,27,22,20],1,'muT')
+distancias=linspace(1,10,10,0.01,'cm') 
+unidade_errada=linspace(1,10,10,0.001,'kg')
+
+def test_regressao_linear_unidades():
+    linha=regressao_linear(distancias,campo_magnético) 
+    linha(distancias)
+    unidade_errada=linspace(1,10,10,0.001,'')
+    with pytest.raises(ValueError):
+       linha(unidade_errada)
+    comparar_medidas(linha(distancias)[0],campo_magnético[0])
+
+def test_regressao_cubica_unidades():
+    for grau in [1,2,3,4,5,6]:
+        cubica=regressao_polinomial(distancias,campo_magnético,grau) 
+        cubica(distancias)
+        with pytest.raises(ValueError):
+            cubica(unidade_errada)
+        comparar_medidas(cubica(distancias)[0],campo_magnético[0])
+
+
+def test_regressao_exponencial_unidades():
+    exponencial=regressao_exponencial(distancias,campo_magnético) 
+    exponencial(distancias)
+    with pytest.raises(pint.errors.DimensionalityError):
+       exponencial(unidade_errada)
+    comparar_medidas(exponencial(distancias)[0],campo_magnético[0])
+
+def test_regressao_potencia_unidades():
+    potencia=regressao_potencia(distancias,campo_magnético) 
+    potencia(distancias)
+    with pytest.raises(ValueError):
+       potencia(unidade_errada)
+    comparar_medidas(potencia(distancias)[0],campo_magnético[0])

tests/test_regressoes_polinomiais.py

@@ -21,10 +21,10 @@ def test_regressao_linear():
         y_dados=lab2.arrayM(y,0.01,'m')
         reta=lab2.regressao_linear(x_dados,y_dados)
         a,b=reta
-        assert np.isclose(a.nominal(""),resposta['a'],rtol=1e-3)
-        assert np.isclose(b.nominal(""),resposta['b'],rtol=1e-3)
-        assert np.isclose(a.incerteza(""),resposta['Δa'],rtol=1e-3)
-        assert np.isclose(b.incerteza(""),resposta['Δb'],rtol=1e-3)
+        assert np.isclose(a.nominal("m/s"),resposta['a'],rtol=1e-3)
+        assert np.isclose(b.nominal("m"),resposta['b'],rtol=1e-3)
+        assert np.isclose(a.incerteza("m/s"),resposta['Δa'],rtol=1e-3)
+        assert np.isclose(b.incerteza("m"),resposta['Δb'],rtol=1e-3)
 def test_regressao_polinominal_nominal():
     num=10
     a,b=np.random.normal(0,1,2)
@@ -36,9 +36,9 @@ def test_regressao_polinominal_nominal():
     y_dados=lab2.arrayM(y,0.01,'m')
     parabola=lab2.regressao_polinomial(x_dados,y_dados,2)
     parabola_predito,a_predito,b_predito=parabola
-    assert np.isclose(parabola_predito.nominal(""),0,atol=1e-2)
-    assert np.isclose(a_predito.nominal(""),a,rtol=1e-1)
-    assert np.isclose(b_predito.nominal(""),b,rtol=1e-1)
+    assert np.isclose(parabola_predito.nominal("m/s²"),0,atol=1e-2)
+    assert np.isclose(a_predito.nominal("m/s"),a,rtol=1e-1)
+    assert np.isclose(b_predito.nominal("m"),b,rtol=1e-1)
 
 def test_regressao_polinomial_basic():
     x_dados = lab2.arrayM([1, 2, 3, 4, 5],0,'')