2012-02-24 8 views
6

Ho di nuovo un problema con il pacchetto vector-space. Ho ricevuto una risposta molto utile da @mnish in un recente post, ma lì ho affrontato solo una funzione che dipende solo da 1 variabile. cosa succede quando ho, per esempio, una funzione che mappa da coordinate polari a cartesianiDerivati ​​di una funzione multivariata e corrispondente Jacobiano con pacchetto spazio vettoriale

f:(0,oo) x [0,2pi] -> R² 
(r,phi) -> (r*cos(phi),r*sin(phi)) 

che dipende da 2 variabili.

Ho provato questo fuori, con un approccio molto ingenuo:

polar :: Double -> Double -> ((Double,Double) :~> (Double,Double)) 
polar r phi = \(r,phi) -> (((idD) r)*cos(idD phi),((idD) r)*sin(idD phi)) 

ottengo il seguente errore:

Couldn't match expected type `(Double, Double) :> (Double, Double)' 
      with actual type `(t0, t1)' 
In the expression: 
    (((idD) r) * cos (idD phi), ((idD) r) * sin (idD phi)) 
In the expression: 
    \ (r, phi) 
    -> (((idD) r) * cos (idD phi), ((idD) r) * sin (idD phi)) 
In an equation for `polar': 
    polar r phi 
     = \ (r, phi) 
      -> (((idD) r) * cos (idD phi), ((idD) r) * sin (idD phi)) 

Per un componente

polarx :: Double -> Double -> ((Double,Double) :~> Double) 
polarx r phi = \(r,phi) -> ((idD) r)*cos(idD phi) 

ottengo

Couldn't match expected type `Double' 
      with actual type `(Double, Double)' 
Expected type: (Double, Double) :> Double 
    Actual type: (Double, Double) :> (Double, Double) 
In the return type of a call of `idD' 
In the first argument of `(*)', namely `((idD) r)' 

Apparentemente c'è qualche tipo di disturbo, ma non riesco a capire cosa c'è che non va.

Un'altra domanda sorge quando voglio calcolare lo Jacobiano di tale mappatura. Come suggerisce il nome, ha qualcosa a che fare con le mappe lineari, che è, ovviamente, coperto dal pacchetto, in realtà si basa su quelle mappe. Ma ancora una volta, la mia conoscenza di Haskell è insufficiente, per ricavare una soluzione per conto mio.

+1

Mi sembra di ricordare che un limite chiave dell'elegantissima formulazione di differenziazione automatica di Conal è che funziona solo su derivati ​​lungo un singolo asse. Se vuoi Jacobians, ecc., Penso che il pacchetto pubblicitario di ekmett sia la strada da percorrere: http://hackage.haskell.org/package/ad-1.3.0.1 – sclv

+1

Grazie a @sclv, ho appena esaminato questo modulo e io devo dire, wow, sono impressionato. Non ho notato questo pacchetto e farò un tentativo, grazie per aver segnalato questo – TheMADMAN

+0

Non sei solo - Sto facendo fatica a capire come i vari tipi multidimensionali combaciano. Vado a leggere il documento "Beautiful Differentiation" e spero che faccia luce: il pacchetto pubblicitario sembra un po 'più semplice sui tipi! – Oliver

risposta

2

Ho finalmente trovato una soluzione al mio problema, non è stato difficile, ma mi ci è voluto un po 'per capirlo. Nel caso in cui qualcun altro sia interessato, vi presento i dettagli.

Primo Ecco il mio codice per il caso polare:

polarCoordD :: ((Double,Double) :~> (Double,Double)) 
polarCoordD = \(r,phi) -> pairD (polarx (r,phi), polary (r,phi)) 
where polarx :: (Double,Double) :~> Double 
     polarx = \(r,phi) -> (fst . unpairD $ (idD) (r,phi))*cos(snd . unpairD $ idD (r, phi)) 
     polary :: (Double,Double) :~> Double 
     polary = \(r,phi) -> (fst . unpairD $ (idD) (r,phi))*sin(snd . unpairD $ idD (r, phi)) 

La chiave è stato quello di rendere la "variabile di derivazione" (idD) consapevoli della tupla (r, phi) che detiene le due variabili che voglio per differenziare. Quindi devo decomprimere la tupla tramite unpairD e scegliere la prima e la seconda parte della coppia risultante (in polarx e polary). Entrambi sono imballati di nuovo in una coppia. Forse c'è un modo più elegante per farlo, ma è così che l'ho capito alla fine.

Da qui non è difficile andare oltre alle coordinate cilindriche o, in effetti, a qualsiasi altro sistema di coordinate ortogonali curve. Per coordinate cilindriche ottengo:

cylCoordD :: (Vec3 Double :~> Vec3 Double) 
cylCoordD = \(r,phi,z) -> tripleD (cylx (r,phi,z), cyly (r,phi,z),cylz (0,0,z)) 
where cylx :: (Double,Double,Double) :~> Double 
     cylx = \(r,phi,z) -> (fst' . untripleD $ (idD) (r,phi,z))*cos(snd' . untripleD $ idD (r, phi,z)) 
     cyly :: (Double,Double,Double) :~> Double 
     cyly = \(r,phi,z) -> (fst' . untripleD $ (idD) (r,phi,z))*sin(snd' . untripleD $ idD (r, phi,z)) 
     cylz :: (Double,Double,Double) :~> Double 
     cylz = \(_,_,z) -> third . untripleD $ idD (0,0,z) 
     fst' :: (a,b,c) -> a 
     fst' (x,_,_) = x 
     snd' :: (a,b,c) -> b 
     snd' (_,y,_) = y 
     third :: (a,b,c) -> c 
     third (_,_,z) = z 

dove Vec3 Double appartiene type Vec3 a = (a, a, a). Ora ci può anche costruire una matrice di trasformazione:

let transmat = \(r,phi,z) -> powVal $ liftD3 (,,) (normalized $ derivAtBasis (cylCoordD (r,phi,z)) (Left())) (normalized $ derivAtBasis (cylCoordD (r,phi,z)) (Right (Left()))) (normalized $ derivAtBasis (cylCoordD (r,phi,z)) (Right (Right()))) 

*Main> transmat (2, rad 0, 0) 
((1.0,0.0,0.0),(0.0,1.0,0.0),(0.0,0.0,1.0)) 

*Main> transmat (2, rad 90, 0) 
((6.123233995736766e-17,1.0,0.0),(-1.0,6.123233995736766e-17,0.0),(0.0,0.0,1.0)) 

rad è una funzione di convenienza

rad :: Double -> Double 
rad = (pi*) . (*recip 180) 

Ora sarebbe interessante per convertire questa "matrice" per il tipo di matrice di Numeric Prelude e/o hmatrix , ma non sono sicuro che sarebbe utile. Tuttavia, sarebbe un buon esempio per l'uso del pacchetto vector-space.

Devo ancora capire l'utilizzo e soprattutto l'applicazione delle mappe lineari.

0

Ho appena visto questa domanda successiva. Non sono sicuro di ciò che si vuole:

  • la matrice Jacobiana
  • un Jacobiano-vettore prodotto
  • un prodotto jacobiano-trasposizione-vettore

In un tale sistema a bassa dimensionalità, Assumerò il primo. (Gli altri tornano utili soprattutto quando il sistema è sufficientemente dimensionato da non voler memorizzare o calcolare la periosa Jacobiana, ma invece trattarla come una matrice sparsa generalizzata.) In ogni caso:

Prelude> :m + Numeric.AD 
Prelude Numeric.AD> let f [r,phi] = map (r*) [cos phi, sin phi] 
Prelude Numeric.AD> jacobian f [2,3::Float] 
[[-0.9899925,-0.28224],[0.14112,-1.979985]] 
Problemi correlati