text:
{-# OPTIONS_GHC -Wall -Wno-partial-type-signatures #-}

{-# LANGUAGE QuasiQuotes #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE MultiWayIf #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE PartialTypeSignatures #-}

-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Sundials.CVode.ODE
-- Copyright   :  Dominic Steinitz 2018,
--                Novadiscovery 2018
-- License     :  BSD
-- Maintainer  :  Dominic Steinitz
-- Stability   :  provisional
--
-- Solution of ordinary differential equation (ODE) initial value problems.
--
-- <https://computation.llnl.gov/projects/sundials/sundials-software>
--
-- A simple example:
--
-- <<diagrams/brusselator.png#diagram=brusselator&height=400&width=500>>
--
-- @
-- import           Numeric.Sundials.CVode.ODE
-- import           Numeric.LinearAlgebra
--
-- import           Plots as P
-- import qualified Diagrams.Prelude as D
-- import           Diagrams.Backend.Rasterific
--
-- brusselator :: Double -> [Double] -> [Double]
-- brusselator _t x = [ a - (w + 1) * u + v * u * u
--                    , w * u - v * u * u
--                    , (b - w) / eps - w * u
--                    ]
--   where
--     a = 1.0
--     b = 3.5
--     eps = 5.0e-6
--     u = x !! 0
--     v = x !! 1
--     w = x !! 2
--
-- lSaxis :: [[Double]] -> P.Axis B D.V2 Double
-- lSaxis xs = P.r2Axis &~ do
--   let ts = xs!!0
--       us = xs!!1
--       vs = xs!!2
--       ws = xs!!3
--   P.linePlot' $ zip ts us
--   P.linePlot' $ zip ts vs
--   P.linePlot' $ zip ts ws
--
-- main = do
--   let res1 = odeSolve brusselator [1.2, 3.1, 3.0] (fromList [0.0, 0.1 .. 10.0])
--   renderRasterific "diagrams/brusselator.png"
--                    (D.dims2D 500.0 500.0)
--                    (renderAxis $ lSaxis $ [0.0, 0.1 .. 10.0]:(toLists $ tr res1))
-- @
--
-----------------------------------------------------------------------------
module Numeric.Sundials.CVode.ODE ( odeSolve
                                   , odeSolveV
                                   , odeSolveVWith
                                   , odeSolveVWith'
                                   , odeSolveRootVWith'
                                   , odeSolveWithEvents
                                   , ODEMethod(..)
                                   , StepControl(..)
                                   , SolverResult(..)
                                   , SundialsSolution(..)
                                   , EventSpec(..)
                                   , EventInfo(..)
                                   , CrossingDirection(..)
                                   ) where

import qualified Language.C.Inline as C
import qualified Language.C.Inline.Unsafe as CU

import           Data.Monoid ((<>))
import           Data.Maybe (isJust, fromJust)
import           Data.List (genericLength)

import           Foreign.C.Types (CDouble, CInt, CLong)
import           Foreign.Ptr (Ptr)
import           Foreign.Storable (peek, poke)

import qualified Data.Vector.Storable as V

import           Data.Coerce (coerce)
import           System.IO.Unsafe (unsafePerformIO)

import           Numeric.LinearAlgebra.Devel (createVector)

import           Numeric.LinearAlgebra.HMatrix (Vector, Matrix, toList, rows,
                                                cols, toLists, size, reshape,
                                                subVector, subMatrix, toColumns)

import           Numeric.Sundials.Arkode (cV_ADAMS, cV_BDF,
                                          vectorToC, cV_SUCCESS,
                                          SunVector(..))
import qualified Numeric.Sundials.Arkode as T
import           Numeric.Sundials.ODEOpts

C.context (C.baseCtx <> C.vecCtx <> C.funCtx <> T.sunCtx)

C.include "<stdlib.h>"
C.include "<stdio.h>"
C.include "<math.h>"
C.include "<cvode/cvode.h>"               -- prototypes for CVODE fcts., consts.
C.include "<nvector/nvector_serial.h>"    -- serial N_Vector types, fcts., macros
C.include "<sunmatrix/sunmatrix_dense.h>" -- access to dense SUNMatrix
C.include "<sunlinsol/sunlinsol_dense.h>" -- access to dense SUNLinearSolver
C.include "<cvode/cvode_direct.h>"        -- access to CVDls interface
C.include "<cvode/cvode_diag.h>"          -- CVDiag
C.include "<sundials/sundials_types.h>"   -- definition of type realtype
C.include "<sundials/sundials_math.h>"
C.include "../../../helpers.h"
C.include "Numeric/Sundials/Arkode_hsc.h"

-- | The direction in which a function should cross the x axis
data CrossingDirection = Upwards | Downwards | AnyDirection
  deriving (Eq, Show)

-- Contrary to the documentation, it appears that CVodeGetRootInfo
-- may use both 1 and -1 to indicate a root, depending on the
-- direction of the sign change. See near the end of cvRootfind.
intToDirection :: Integral d => d -> Maybe CrossingDirection
intToDirection d =
  case d of
    1  -> Just Upwards
    -1 -> Just Downwards
    _  -> Nothing

-- | Almost inverse of 'intToDirection'. Map 'Upwards' to 1, 'Downwards' to
-- -1, and 'AnyDirection' to 0.
directionToInt :: Integral d => CrossingDirection -> d
directionToInt d =
  case d of
    Upwards -> 1
    Downwards -> -1
    AnyDirection -> 0

data SundialsSolution =
  SundialsSolution
  { actualTimeGrid :: V.Vector Double                 -- ^ actual time grid returned by the solver (with duplicated event times)
  , solutionMatrix :: Matrix Double                -- ^ matrix of solutions: each column is an unknwown (add also the time vector?)
  , eventInfo      :: [EventInfo]         -- ^ event infos, as many items as triggered events during the simulation
  , diagnostics    :: SundialsDiagnostics -- ^ usual Sundials diagnostics
  }

data EventInfo =
  EventInfo
  { eventTime     :: !Double            -- ^ time at which event was triggered
  , eventIndex    :: !Int               -- ^ which index was triggered
  , rootDirection :: !CrossingDirection -- ^ in which direction ((+)->(-) or (-)->(+)) the root is crossed
  }
  deriving Show

getMethod :: ODEMethod -> Int
getMethod (ADAMS) = cV_ADAMS
getMethod (BDF)   = cV_BDF

getJacobian :: ODEMethod -> Maybe Jacobian
getJacobian _ = Nothing

-- | A version of 'odeSolveVWith' with reasonable default step control.
odeSolveV
    :: ODEMethod
    -> Maybe Double      -- ^ initial step size - by default, CVode
                         -- estimates the initial step size to be the
                         -- solution \(h\) of the equation
                         -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
                         -- \(\ddot{y}\) is an estimated value of the
                         -- second derivative of the solution at \(t_0\)
    -> Double            -- ^ absolute tolerance for the state vector
    -> Double            -- ^ relative tolerance for the state vector
    -> (Double -> Vector Double -> Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
    -> Vector Double     -- ^ initial conditions
    -> Vector Double     -- ^ desired solution times
    -> Matrix Double     -- ^ solution
odeSolveV meth hi epsAbs epsRel f y0 ts =
  odeSolveVWith meth (X epsAbs epsRel) hi g y0 ts
  where
    g t x0 = coerce $ f t x0

-- | A version of 'odeSolveV' with reasonable default parameters and
-- system of equations defined using lists. FIXME: we should say
-- something about the fact we could use the Jacobian but don't for
-- compatibility with hmatrix-gsl.
odeSolve :: (Double -> [Double] -> [Double]) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
         -> [Double]                         -- ^ initial conditions
         -> Vector Double                    -- ^ desired solution times
         -> Matrix Double                    -- ^ solution
odeSolve f y0 ts =
  -- FIXME: These tolerances are different from the ones in GSL
  odeSolveVWith BDF (XX' 1.0e-6 1.0e-10 1 1)  Nothing g (V.fromList y0) (V.fromList $ toList ts)
  where
    g t x0 = V.fromList $ f t (V.toList x0)

-- | A version of 'odeSolveVWith'' with reasonable default solver
-- options.
odeSolveVWith ::
  ODEMethod
  -> StepControl
  -> Maybe Double -- ^ initial step size - by default, CVode
                  -- estimates the initial step size to be the
                  -- solution \(h\) of the equation
                  -- \(\|\frac{h^2\ddot{y}}{2}\| = 1\), where
                  -- \(\ddot{y}\) is an estimated value of the second
                  -- derivative of the solution at \(t_0\)
  -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
  -> V.Vector Double                     -- ^ Initial conditions
  -> V.Vector Double                     -- ^ Desired solution times
  -> Matrix Double                       -- ^ Error code or solution
odeSolveVWith method control initStepSize f y0 tt =
  case odeSolveVWith' opts f y0 tt of
    Left  (c, _v) -> error $ show c -- FIXME
    Right (v, _d) -> v
  where
    opts = ODEOpts { maxNumSteps = 10000
                   , minStep     = 1.0e-12
                   , maxFail     = 10
                   , odeMethod   = method
                   , stepControl = control
                   , initStep    = initStepSize
                   }

odeSolveVWith' ::
  ODEOpts
  -> (Double -> V.Vector Double -> V.Vector Double) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
  -> V.Vector Double                     -- ^ Initial conditions
  -> V.Vector Double                     -- ^ Desired solution times
  -> Either (Matrix Double, Int) (Matrix Double, SundialsDiagnostics) -- ^ Error code or solution
odeSolveVWith' opts f y0 tt =
  case solveOdeC (fromIntegral $ maxFail opts)
                  (fromIntegral $ maxNumSteps opts) (coerce $ minStep opts)
                  (fromIntegral . getMethod . odeMethod $ opts) (coerce $ initStep opts) jacH (scise $ stepControl opts)
                  (coerce f) (coerce y0)
                  0 (\_ x -> x) [] 0 (\_ _ y -> y) (coerce tt) of
    -- Remove the time column for backwards compatibility
    SolverError m c         -> Left
                               ( subMatrix (0, 1) (V.length tt, l) m
                               , fromIntegral c
                               )
    SolverSuccess _ m d     -> Right
                               ( subMatrix (0, 1) (V.length tt, l) m
                               , d
                               )
  where
    l = size y0
    scise (X aTol rTol)                          = coerce (V.replicate l aTol, rTol)
    scise (X' aTol rTol)                         = coerce (V.replicate l aTol, rTol)
    scise (XX' aTol rTol yScale _yDotScale)      = coerce (V.replicate l aTol, yScale * rTol)
    -- FIXME; Should we check that the length of ss is correct?
    scise (ScXX' aTol rTol yScale _yDotScale ss) = coerce (V.map (* aTol) ss, yScale * rTol)
    jacH = fmap (\g t v -> matrixToSunMatrix $ g (coerce t) (coerce v)) $
           getJacobian $ odeMethod opts

matrixToSunMatrix :: Matrix Double -> T.SunMatrix
matrixToSunMatrix m = T.SunMatrix { T.rows = nr, T.cols = nc, T.vals = vs }
  where
    nr = fromIntegral $ rows m
    nc = fromIntegral $ cols m
    -- FIXME: efficiency
    vs = V.fromList $ map coerce $ concat $ toLists m

solveOdeC ::
  CInt ->
  CLong ->
  CDouble ->
  CInt ->
  Maybe CDouble ->
  (Maybe (CDouble -> V.Vector CDouble -> T.SunMatrix)) ->
  (V.Vector CDouble, CDouble) ->
  (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
  -> V.Vector CDouble -- ^ Initial conditions
  -> CInt -- ^ Number of event equations
  -> (CDouble -> V.Vector CDouble -> V.Vector CDouble) -- ^ The event equations themselves
  -> [CrossingDirection] -- ^ The required crossing direction for each event
  -> CInt -- ^ Maximum number of events
  -> (Int -> CDouble -> V.Vector CDouble -> V.Vector CDouble)
      -- ^ Function to reset/update the state when an event occurs. The
      -- 'Int' argument is the 0-based number of the event that has
      -- occurred. If multiple events have occurred at the same time, they
      -- are handled in the increasing order of the event index. The other
      -- arguments are the time and the point in the state space. Return
      -- the updated point in the state space.
  -> V.Vector CDouble -- ^ Desired solution times
  -> SolverResult
solveOdeC maxErrTestFails maxNumSteps_ minStep_ method initStepSize
          jacH (aTols, rTol) fun f0 nr event_fn directions max_events apply_event ts =
  unsafePerformIO $ do

  let isInitStepSize :: CInt
      isInitStepSize = fromIntegral $ fromEnum $ isJust initStepSize
      ss :: CDouble
      ss = case initStepSize of
             -- It would be better to put an error message here but
             -- inline-c seems to evaluate this even if it is never
             -- used :(
             Nothing -> 0.0
             Just x  -> x

  let dim = V.length f0
      nEq :: CLong
      nEq = fromIntegral dim
      nTs :: CInt
      nTs = fromIntegral $ V.length ts
  output_mat_mut :: V.MVector _ CDouble <- V.thaw =<< createVector ((1 + fromIntegral dim) * (fromIntegral (2 * max_events) + fromIntegral nTs))
  diagMut :: V.MVector _ CLong <- V.thaw =<< createVector 10 -- FIXME
  -- We need the types that sundials expects.
  -- FIXME: The Haskell type is currently empty!
  let funIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr () -> IO CInt
      funIO t y f _ptr = do
        sv <- peek y
        poke f $ SunVector { sunVecN = sunVecN sv
                           , sunVecVals = fun t (sunVecVals sv)
                           }
        return 0

  let nrPre = fromIntegral nr
  gResults :: V.Vector CInt <- createVector nrPre
  -- FIXME: Do we need to do this here? Maybe as it will get GC'd and
  -- we'd have to do a malloc in C otherwise :(
  gResMut <- V.thaw gResults
  event_index_mut :: V.MVector _ CInt <- V.thaw =<< createVector (fromIntegral max_events)
  event_time_mut :: V.MVector _ CDouble <- V.thaw =<< createVector (fromIntegral max_events)
  -- Total number of events. This is *not* directly re
  n_events_mut :: V.MVector _ CInt <- V.thaw =<< createVector 1
  -- Total number of rows in the output_mat_mut matrix. It *cannot* be
  -- inferred from n_events_mut because when an event occurs k times, it
  -- contributes k to n_events_mut but only 2 to n_rows_mut.
  n_rows_mut :: V.MVector _ CInt <- V.thaw =<< createVector 1
  actual_event_direction_mut :: V.MVector _ CInt <- V.thaw =<< createVector (fromIntegral max_events)

  let event_fn_c :: CDouble -> Ptr T.SunVector -> Ptr CDouble -> Ptr () -> IO CInt
      event_fn_c x y f _ptr = do
        vals <- event_fn x <$> (sunVecVals <$> peek y)
        -- FIXME: We should be able to use poke somehow
        vectorToC vals nrPre f
        return 0

      apply_event_c :: CInt -> CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> IO CInt
      apply_event_c event_index t y y' = do
        sv <- peek y
        poke y' $ SunVector
          { sunVecN = sunVecN sv
          , sunVecVals = apply_event (fromIntegral event_index) t (sunVecVals sv)
          }
        return 0

      requested_event_directions :: V.Vector CInt
      requested_event_directions = V.fromList $ map directionToInt directions

  let isJac :: CInt
      isJac = fromIntegral $ fromEnum $ isJust jacH
      jacIO :: CDouble -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunMatrix ->
               Ptr () -> Ptr T.SunVector -> Ptr T.SunVector -> Ptr T.SunVector ->
               IO CInt
      jacIO t y _fy jacS _ptr _tmp1 _tmp2 _tmp3 = do
        case jacH of
          Nothing   -> error "Numeric.Sundials.CVode.ODE: Jacobian not defined"
          Just jacI -> do j <- jacI t <$> (sunVecVals <$> peek y)
                          poke jacS j
                          -- FIXME: I don't understand what this comment means
                          -- Unsafe since the function will be called many times.
                          [CU.exp| int{ 0 } |]

  res :: Int <- fromIntegral <$> [C.block| int {
                         /* general problem variables */

                         int flag;                  /* reusable error-checking flag                 */

                         int i, j;                  /* reusable loop indices                        */
                         N_Vector y = NULL;         /* empty vector for storing solution            */
                         N_Vector tv = NULL;        /* empty vector for storing absolute tolerances */

                         SUNMatrix A = NULL;        /* empty matrix for linear solver               */
                         SUNLinearSolver LS = NULL; /* empty linear solver object                   */
                         void *cvode_mem = NULL;    /* empty CVODE memory structure                 */
                         realtype t;
                         long int nst, nfe, nsetups, nje, nfeLS, nni, ncfn, netf, nge;

                         realtype tout;

                         /* input_ind tracks the current index into the ts array */
                         int input_ind = 1;
                         /* output_ind tracks the current row into the output_mat_mut matrix.
                            If differs from input_ind because of the extra rows corresponding to events. */
                         int output_ind = 1;
                         /* We need to update n_rows_mut every time we update output_ind because
                            of the possibility of early return (in which case we still need to assemble
                            the partial results matrix). We could even work with n_rows_mut only and ditch
                            output_ind, but the inline-c expression is quite verbose, and output_ind is
                            more convenient to use in index calculations.
                         */
			 ($vec-ptr:(int *n_rows_mut))[0] = output_ind;
                         /* event_ind tracks the current event number */
                         int event_ind = 0;

                         /* general problem parameters */

                         realtype T0 = RCONST(($vec-ptr:(double *ts))[0]); /* initial time              */
                         sunindextype NEQ = $(sunindextype nEq);           /* number of dependent vars. */

                         /* Initialize data structures */

                         y = N_VNew_Serial(NEQ); /* Create serial vector for solution */
                         if (check_flag((void *)y, "N_VNew_Serial", 0)) return 1;
                         /* Specify initial condition */
                         for (i = 0; i < NEQ; i++) {
                           NV_Ith_S(y,i) = ($vec-ptr:(double *f0))[i];
                         };

                         cvode_mem = CVodeCreate($(int method), CV_NEWTON);
                         if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);

                         /* Call CVodeInit to initialize the integrator memory and specify the
                          * user's right hand side function in y'=f(t,y), the inital time T0, and
                          * the initial dependent variable vector y. */
                         flag = CVodeInit(cvode_mem, $fun:(int (* funIO) (double t, SunVector y[], SunVector dydt[], void * params)), T0, y);
                         if (check_flag(&flag, "CVodeInit", 1)) return(1);

                         tv = N_VNew_Serial(NEQ); /* Create serial vector for absolute tolerances */
                         if (check_flag((void *)tv, "N_VNew_Serial", 0)) return 1;
                         /* Specify tolerances */
                         for (i = 0; i < NEQ; i++) {
                           NV_Ith_S(tv,i) = ($vec-ptr:(double *aTols))[i];
                         };

                         flag = CVodeSetMinStep(cvode_mem, $(double minStep_));
                         if (check_flag(&flag, "CVodeSetMinStep", 1)) return 1;
                         flag = CVodeSetMaxNumSteps(cvode_mem, $(long int maxNumSteps_));
                         if (check_flag(&flag, "CVodeSetMaxNumSteps", 1)) return 1;
                         flag = CVodeSetMaxErrTestFails(cvode_mem, $(int maxErrTestFails));
                         if (check_flag(&flag, "CVodeSetMaxErrTestFails", 1)) return 1;

                         /* Call CVodeSVtolerances to specify the scalar relative tolerance
                          * and vector absolute tolerances */
                         flag = CVodeSVtolerances(cvode_mem, $(double rTol), tv);
                         if (check_flag(&flag, "CVodeSVtolerances", 1)) return(1);

                         /* Call CVodeRootInit to specify the root function event_fn_c with nr components */
                         flag = CVodeRootInit(cvode_mem, $(int nr), $fun:(int (* event_fn_c) (double t, SunVector y[], double gout[], void * params)));

                         if (check_flag(&flag, "CVodeRootInit", 1)) return(1);

                         /* Set the initial step size if there is one */
                         if ($(int isInitStepSize)) {
                           /* FIXME: We could check if the initial step size is 0 */
                           /* or even NaN and then throw an error                 */
                           flag = CVodeSetInitStep(cvode_mem, $(double ss));
                           if (check_flag(&flag, "CVodeSetInitStep", 1)) return 1;
                         }

                         /* Set the Jacobian if there is one */
                         if ($(int isJac)) {
                           /* Initialize dense matrix data structure and solver */
                           A = SUNDenseMatrix(NEQ, NEQ);
                           if (check_flag((void *)A, "SUNDenseMatrix", 0)) return 1;
                           LS = SUNDenseLinearSolver(y, A);
                           if (check_flag((void *)LS, "SUNDenseLinearSolver", 0)) return 1;

                           /* Attach matrix and linear solver */
                           flag = CVDlsSetLinearSolver(cvode_mem, LS, A);
                           if (check_flag(&flag, "CVDlsSetLinearSolver", 1)) return 1;

                           /* Set the Jacobian */
                           flag = CVDlsSetJacFn(cvode_mem, $fun:(int (* jacIO) (double t, SunVector y[], SunVector fy[], SunMatrix Jac[], void * params, SunVector tmp1[], SunVector tmp2[], SunVector tmp3[])));
                           if (check_flag(&flag, "CVDlsSetJacFn", 1)) return 1;
                         } else /* Else we use CVDiag directly */ {
                           flag = CVDiag(cvode_mem);
                           if (check_flag(&flag, "CVDiag", 1)) return 1;
                         }

                         /* Store initial conditions */
			 ($vec-ptr:(double *output_mat_mut))[0 * (NEQ + 1) + 0] = ($vec-ptr:(double *ts))[0];
                         for (j = 0; j < NEQ; j++) {
                           ($vec-ptr:(double *output_mat_mut))[0 * (NEQ + 1) + (j + 1)] = NV_Ith_S(y,j);
                         }

                         while (1) {
                           flag = CVode(cvode_mem, ($vec-ptr:(double *ts))[input_ind], y, &t, CV_NORMAL); /* call integrator */
                           if (check_flag(&flag, "CVode solver failure, stopping integration", 1)) return 1;

                           /* Store the results for Haskell */
			   ($vec-ptr:(double *output_mat_mut))[output_ind * (NEQ + 1) + 0] = t;
                           for (j = 0; j < NEQ; j++) {
                             ($vec-ptr:(double *output_mat_mut))[output_ind * (NEQ + 1) + (j + 1)] = NV_Ith_S(y,j);
                           }
                           output_ind++;
                           ($vec-ptr:(int *n_rows_mut))[0] = output_ind;

                           if (flag == CV_ROOT_RETURN) {
                             if (event_ind >= $(int max_events)) {
                               /* We don't have any more space for events. Return an error. */
                               return 1;
                             }

                             /* Are we interested in this event?
                                If not, continue without any observable side-effects.
                             */
                             int good_event = 0;
			     flag = CVodeGetRootInfo(cvode_mem, ($vec-ptr:(int *gResMut)));
			     if (check_flag(&flag, "CVodeGetRootInfo", 1)) return 1;
                             for (i = 0; i < $(int nr); i++) {
                               int ev = ($vec-ptr:(int *gResMut))[i];
                               int req_dir = ($vec-ptr:(const int *requested_event_directions))[i];
                               if (ev != 0 && ev * req_dir >= 0) {
                                 good_event = 1;

                                 ($vec-ptr:(int *actual_event_direction_mut))[event_ind] = ev;
                                 ($vec-ptr:(int *event_index_mut))[event_ind] = i;
                                 ($vec-ptr:(double *event_time_mut))[event_ind] = t;
                                 event_ind++;

                                 /* Update the state with the supplied function */
                                 $fun:(int (* apply_event_c) (int, double, SunVector y[], SunVector z[]))(i, t, y, y);
                               }
                             }

                             if (good_event) {
                               ($vec-ptr:(double *output_mat_mut))[output_ind * (NEQ + 1) + 0] = t;
                               for (j = 0; j < NEQ; j++) {
                                 ($vec-ptr:(double *output_mat_mut))[output_ind * (NEQ + 1) + (j + 1)] = NV_Ith_S(y,j);
                               }
                               output_ind++;
                               ($vec-ptr:(int *n_rows_mut))[0] = output_ind;

                               flag = CVodeReInit(cvode_mem, t, y);
                               if (check_flag(&flag, "CVodeReInit", 1)) return(1);
                             } else {
                               /* Since this is not a wanted event, it shouldn't get a row */
                               output_ind--;
                               ($vec-ptr:(int *n_rows_mut))[0] = output_ind;
                             }
                           }
			   else {
			     if (++input_ind >= $(int nTs))
                               break;
			   }
                         }

			 /* The number of actual roots we found */
			 ($vec-ptr:(int *n_events_mut))[0] = event_ind;

                         /* Get some final statistics on how the solve progressed */
                         flag = CVodeGetNumSteps(cvode_mem, &nst);
                         check_flag(&flag, "CVodeGetNumSteps", 1);
                         ($vec-ptr:(long int *diagMut))[0] = nst;

                         /* FIXME */
                         ($vec-ptr:(long int *diagMut))[1] = 0;

                         flag = CVodeGetNumRhsEvals(cvode_mem, &nfe);
                         check_flag(&flag, "CVodeGetNumRhsEvals", 1);
                         ($vec-ptr:(long int *diagMut))[2] = nfe;
                         /* FIXME */
                         ($vec-ptr:(long int *diagMut))[3] = 0;

                         flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);
                         check_flag(&flag, "CVodeGetNumLinSolvSetups", 1);
                         ($vec-ptr:(long int *diagMut))[4] = nsetups;

                         flag = CVodeGetNumErrTestFails(cvode_mem, &netf);
                         check_flag(&flag, "CVodeGetNumErrTestFails", 1);
                         ($vec-ptr:(long int *diagMut))[5] = netf;

                         flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);
                         check_flag(&flag, "CVodeGetNumNonlinSolvIters", 1);
                         ($vec-ptr:(long int *diagMut))[6] = nni;

                         flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);
                         check_flag(&flag, "CVodeGetNumNonlinSolvConvFails", 1);
                         ($vec-ptr:(long int *diagMut))[7] = ncfn;

                         flag = CVDlsGetNumJacEvals(cvode_mem, &nje);
                         check_flag(&flag, "CVDlsGetNumJacEvals", 1);
                         ($vec-ptr:(long int *diagMut))[8] = ncfn;

                         flag = CVDlsGetNumRhsEvals(cvode_mem, &nfeLS);
                         check_flag(&flag, "CVDlsGetNumRhsEvals", 1);
                         ($vec-ptr:(long int *diagMut))[9] = ncfn;

                         /* Clean up and return */

                         N_VDestroy(y);          /* Free y vector          */
                         N_VDestroy(tv);         /* Free tv vector         */
                         CVodeFree(&cvode_mem);  /* Free integrator memory */
                         SUNLinSolFree(LS);      /* Free linear solver     */
                         SUNMatDestroy(A);       /* Free A matrix          */

                         return CV_SUCCESS;
                       } |]
  preD <- V.freeze diagMut
  let d = SundialsDiagnostics (fromIntegral $ preD V.!0)
                              (fromIntegral $ preD V.!1)
                              (fromIntegral $ preD V.!2)
                              (fromIntegral $ preD V.!3)
                              (fromIntegral $ preD V.!4)
                              (fromIntegral $ preD V.!5)
                              (fromIntegral $ preD V.!6)
                              (fromIntegral $ preD V.!7)
                              (fromIntegral $ preD V.!8)
                              (fromIntegral $ preD V.!9)
  n_rows <- fromIntegral . V.head <$> V.freeze n_rows_mut
  output_mat <- coerce . reshape (dim + 1) . subVector 0 ((dim + 1) * n_rows) <$>
    V.freeze output_mat_mut
  n_events <- fromIntegral . V.head <$> V.freeze n_events_mut
  event_time             :: V.Vector Double
    <- coerce . V.take n_events <$> V.freeze event_time_mut
  event_index            :: V.Vector Int
    <- V.map fromIntegral . V.take n_events <$> V.freeze event_index_mut
  actual_event_direction :: V.Vector CInt
    <- V.take n_events <$> V.freeze actual_event_direction_mut
  let
    events :: [EventInfo]
    events = zipWith3 EventInfo
      (V.toList event_time)
      (V.toList event_index)
      (map (fromJust . intToDirection) $ V.toList actual_event_direction)
  return $
    if res == cV_SUCCESS
      then
        SolverSuccess events output_mat d
      else
        SolverError output_mat res

data SolverResult
  = SolverError !(Matrix Double) !Int
      -- ^ Partial results and error code
  | SolverSuccess
      [EventInfo]
      !(Matrix Double)
      !SundialsDiagnostics
      -- ^ Times at which the event was triggered, information about which root and the
                                                   -- results and diagnostics.
    deriving Show

odeSolveRootVWith' ::
  ODEOpts
  -> (Double -> V.Vector Double -> V.Vector Double)
      -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
  -> Maybe (Double -> Vector Double -> Matrix Double)
      -- ^ The Jacobian (optional)
  -> V.Vector Double                      -- ^ Initial conditions
  -> [EventSpec]                          -- ^ Event specifications
  -> Int                                  -- ^ Maximum number of events
  -> V.Vector Double                      -- ^ Desired solution times
  -> SolverResult
odeSolveRootVWith' opts f mb_jacobian y0 event_specs nRootEvs tt =
  solveOdeC (fromIntegral $ maxFail opts)
                 (fromIntegral $ maxNumSteps opts) (coerce $ minStep opts)
                 (fromIntegral . getMethod . odeMethod $ opts) (coerce $ initStep opts) jacH (scise $ stepControl opts)
                 (coerce f) (coerce y0)
                 (genericLength event_specs) event_equations event_directions
                 (fromIntegral nRootEvs) reset_state
                 (coerce tt)
  where
    l = size y0
    scise (X aTol rTol)                          = coerce (V.replicate l aTol, rTol)
    scise (X' aTol rTol)                         = coerce (V.replicate l aTol, rTol)
    scise (XX' aTol rTol yScale _yDotScale)      = coerce (V.replicate l aTol, yScale * rTol)
    -- FIXME; Should we check that the length of ss is correct?
    scise (ScXX' aTol rTol yScale _yDotScale ss) = coerce (V.map (* aTol) ss, yScale * rTol)
    jacH = fmap (\g t v -> matrixToSunMatrix $ g (coerce t) (coerce v)) $ mb_jacobian
    event_equations :: CDouble -> Vector CDouble -> Vector CDouble
    event_equations t y = V.fromList $
      map (\ev -> coerce (eventCondition ev) t y) event_specs
    event_directions :: [CrossingDirection]
    event_directions = map eventDirection event_specs
    reset_state :: Int -> CDouble -> Vector CDouble -> Vector CDouble
    reset_state n_event = coerce $ eventUpdate (event_specs !! n_event)

data EventSpec = EventSpec
  { eventCondition :: Double -> V.Vector Double -> Double
  , eventDirection :: CrossingDirection
  , eventUpdate :: Double -> V.Vector Double -> V.Vector Double
  }

odeSolveWithEvents
  :: ODEOpts
  -> [EventSpec]
    -- ^ Event specifications
  -> Int
    -- ^ Maximum number of events
  -> (Double -> V.Vector Double -> V.Vector Double)
    -- ^ The RHS of the system \(\dot{y} = f(t,y)\)
  -> Maybe (Double -> Vector Double -> Matrix Double)
    -- ^ The Jacobian (optional)
  -> V.Vector Double
    -- ^ Initial conditions
  -> V.Vector Double
    -- ^ Desired solution times
  -> Either Int SundialsSolution
    -- ^ Either an error code or a solution
odeSolveWithEvents opts event_specs max_events rhs mb_jacobian initial sol_times =
  let
    result :: SolverResult
    result =
      odeSolveRootVWith' opts rhs mb_jacobian initial event_specs
        max_events sol_times
  in
    case result of
      SolverError _ code -> Left code
      SolverSuccess events mx diagn ->
        Right $ SundialsSolution
            { actualTimeGrid = extractTimeGrid mx
            , solutionMatrix = mx
            , eventInfo = events
            , diagnostics = diagn
            }
  where
    -- The time grid is the first column of the result matrix
    extractTimeGrid :: Matrix Double -> Vector Double
    extractTimeGrid = head . toColumns

max_stars_repo_name: feuerbach/hmatrix-sundials

max_issues_repo_name: feuerbach/hmatrix-sundials

max_forks_repo_name: feuerbach/hmatrix-sundials


text:
module ComplexPlotter where

import Data.Complex
import System.IO()
import System.Process

-- Converts a number to a complex number
toComplex n = n :+ 0

-- Adds two complex numbers
cadd (a :+ b) (c :+ d) = (a + c) :+ (b + d)

-- Subtracts two complex numbers
csub (a :+ b) (c :+ d) = (a - c) :+ (b - d)

-- Multiplies two complex numbers
cmul (a :+ b) (c :+ d) = (a*c - b*d) :+ (a*d + b*c)

---- Complex numbers
i = 0.0 :+ 1.0
zero = 0.0 :+ 0.0
one = 1.0 :+ 0.0
two = 2.0 :+ 0.0

-- Complex functions
-- x^2
square z = cmul z z

-- Uses partial application to return a function g(x, c) = c + f(x)
cPlus c = cadd (c :+ 0.0)

-- x^2 + 25
squarePlus25 z = cadd (25 :+ 0) $ cmul z z

-- x = i
imagPlane _ = i

-- 2x**2 + 2x + 2
complexPoly z = cadd ((two `cmul` z) `cadd` two) (cmul two (square z))

format (a, b, c, d) = show a ++"\t" ++ show b ++"\t" ++ show c ++ "\t" ++ show d ++ "\n\n"

plot f ss a b = do
          let axis = [a, a+ss..b]
          let graph = [(x, y, magnitude $ f (x :+ y), phase $ f (x :+ y)) | x <- axis, y <- axis]

          writeFile "data.txt" $ foldr (++) "" $ map format graph

          gnuplot <- callCommand "gnuplot \"Surface-Magnitude.plt\""
          return gnuplot

{--
Changes the parameter c of a complex function f(x) + c in the interval [a, b] with a given
stepsize ss.

The method returns a list all functions where c is set according to the interval above.
--}

functionList ss a b = [cPlus c | c <- [a, a+ss..b]]

{--

This method gets a lists of functions and from that generates frames that get encoded into a video

Default parameters:
Function: square
Step Size: 0.125
Interval: [-10, 10]

Change according to your needs in the function below.

Example of how to run the function:

*> animate $ functionList 1 (-25) 25

This will create a file "plot.mp4" showing the complex function z^2 + 0 ... z^2 + 25

--}

animate fx = animate' fx 0 where
  animate' [] n = do
                  video <- callCommand $ "ffmpeg -r 30 -f image2 -s 1920x1080 -start_number 0 -i ./pics/frame%d.png -vframes " ++ show (n-1) ++ " -vcodec libx264 -crf 15  -pix_fmt yuv420p plot.mp4"
                  return video

  animate' (f:fs) n = do
                    current <- plot (f . square) 0.125 (-10) (10)
                    return current

                    rename <- callCommand $ "ren \"pics\\output.png\" \"frame" ++ show n ++ ".png"
                    return rename

                    animate' fs (n+1)

ext: hs


text:
{-# LANGUAGE FlexibleInstances #-}
module SchemeParser.Types where

import Data.IORef
import Data.Complex (Complex)
import Data.Ratio (Rational)

import Control.Monad.Reader
import Control.Monad.Except

import qualified Data.Map as M

import Text.Parsec

-- parser

type LispParser = Parsec String ()

data LispVal = LAtom String
             | LNil
             | LString String
             | LBool Bool
             | LChar Char
             | LNumber Integer
             | LFloat Float
             | LRatio Rational
             | LComplex (Complex Double)
             | LList [LispVal]
             | LDottedList [LispVal] LispVal
             | LVector [LispVal]
             | LPrimFunc ([LispVal] -> Scheme LispVal)
             | LFunc { params :: [String], vararg :: Maybe String,
                       closure :: M.Map String LispVal, body :: [LispVal]}
             deriving Eq

instance Eq ([LispVal] -> Scheme LispVal) where
  _ /= _ = False

-- errors

data LispError = NumArgs Integer [LispVal]
               | TypeMismatch String LispVal
               | Parser ParseError
               | BadSpecialForm String LispVal
               | NotFunction String String
               | UnboundVar String String
               | Default String

type ThrowsError = Either LispError

-- environment

type Env_ = IORef (M.Map String (IORef LispVal))

data Env a = Env
  { globalEnv :: a
  , localEnv  :: a
  }

data EnvType = GlobalEnv | LocalEnv

type SchemeEnv = Env Env_

getEnv :: EnvType -> SchemeEnv -> Env_
getEnv GlobalEnv = globalEnv
getEnv LocalEnv = localEnv

setEnv :: EnvType -> SchemeEnv -> Env_ -> SchemeEnv
setEnv GlobalEnv e e_ = e {globalEnv = e_}
setEnv LocalEnv e e_ = e {localEnv = e_}

pureEnv :: Env_ -> IO (M.Map String LispVal)
pureEnv envR = do
  env <- readIORef envR
  let m = M.toList env
  return $ M.empty

nullEnv :: IO SchemeEnv
nullEnv = Env <$> newIORef M.empty <*> newIORef M.empty

testEnv :: IO SchemeEnv
testEnv = do
  envR <- nullEnv
  let e = globalEnv envR
  env <- readIORef e
  valueRef <- newIORef (LNumber 1)
  writeIORef e (M.insert "foo" valueRef env)
  return envR

-- global type

type Scheme = ReaderT SchemeEnv (ExceptT LispError IO)

runScheme :: Scheme a -> SchemeEnv -> IO (Either LispError a)
runScheme sa env = runExceptT $ runReaderT sa env

testScheme :: Scheme a -> IO (Either LispError a)
testScheme sa = do
  e <- testEnv
  runScheme sa e

max_stars_repo_name: alexpeits/haskell-lab-old

max_issues_repo_name: alexpeits/haskell-lab-old

max_forks_repo_name: alexpeits/haskell-lab-old


text:
{-# OPTIONS_GHC  -fno-warn-unused-binds -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-missing-signatures #-}
{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses, UndecidableInstances, FlexibleContexts, TypeSynonymInstances #-}

---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------
-- |
-- | Module : Bench mark morris
-- | Creator: Xiao Ling
-- | Created: 12/8/2015
-- | see    : https://github.com/snoyberg/conduit
-- | TODO   : move this to benchmark directory
-- |
---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------

module Bmorris where

import Data.List
import Data.Random
import Data.Conduit
import Data.List.Split
import qualified Data.Conduit.List as Cl

import Control.Monad.Identity

import Criterion.Main

import Core
import Statistics
import Morris

{-----------------------------------------------------------------------------
  Benchmark  counter vs counter'

    Benchmark results: it's unclear whether counter is faster than counter'
    due to noise:

    benchmarking count vs count'/count 500
      time                 6.527 s    (-1.716 s .. 19.59 s)
                           0.692 R²   (0.283 R² .. 1.000 R²)
      mean                 4.232 s    (2.136 s .. 5.686 s)
      std dev              2.186 s    (0.0 s .. 2.519 s)
      variance introduced by outliers: 74% (severely inflated)

    benchmarking count vs count'/count' 500
      time                 3.315 s    (-6.152 s .. 16.65 s)
                           0.357 R²   (0.239 R² .. 1.000 R²)
      mean                 5.474 s    (2.826 s .. 6.971 s)
      std dev              2.351 s    (0.0 s .. 2.592 s)
      variance introduced by outliers: 74% (severely inflated)

    The variance actually went up as # of counters increased

------------------------------------------------------------------------------}

main :: IO ()
main = defaultMain . return $ bgroup "flat vs nested morris" [
        bench "flat list e=d=0.05"   . nfIO $ morris 0.05 0.05 [1..10000]
      , bench "nested list e=d=0.05" . nfIO $ morris' 0.05 0.05 [1..10000]
  ]

_count cs  = runRVar (Cl.sourceList [1..1000] $$ count cs) StdRandom
_count' cs = runRVar (Cl.sourceList [1..1000] $$ count' cs) StdRandom

(t1,m1) = (10,3)   :: (Int,Int)
(t2,m2) = (100,3)  :: (Int,Int)
(t3,m3) = (8000,3)  :: (Int,Int)

cs' :: Num a => Int -> Int -> [[Counter]]
cs' t m  = replicate m $ replicate t 0

cs :: Num a => Int -> Int -> [Counter]
cs t m = replicate (t*m) 0

-- * conclude: morris 2x faster than morris', not surprisingly
--main :: IO ()
--main = defaultMain . return $ bgroup "morris vs morris'" [
--        bench "morris 10"  . nfIO $ morris  0.5 0.5 [1..10]
--      , bench "morris' 10" . nfIO $ morris' 0.5 0.5 [1..10]

--      , bench "morris 100"  . nfIO $ morris  0.5 0.5 [1..100]
--      , bench "morris' 100" . nfIO $ morris' 0.5 0.5 [1..100]

--      , bench "morris 1000"  . nfIO $ morris  0.5 0.5 [1..1000]
--      , bench "morris' 1000" . nfIO $ morris' 0.5 0.5 [1..1000]

--      , bench "morris 10000"  . nfIO $ morris  0.5 0.5 [1..10000]
--      , bench "morris' 10000" . nfIO $ morris' 0.5 0.5 [1..10000]
--  ]

{-----------------------------------------------------------------------------
  Benchmark  t of m vs m of t
------------------------------------------------------------------------------}

---- * conclusion : counters faster by 1 ms on 0.05 0.05
--main :: IO ()
--main = defaultMain . return $ bgroup "t of m vs m of t tests" [
--        bench "eps = 0.5 d = 0.5 counter " . nfIO $ incrsIO $ counters  0.5 0.5
--      , bench "eps = 0.5 d = 0.5 counter'" . nfIO $ incrsIO $ counters' 0.5 0.5

--      , bench "eps = 0.1 d = 0.1 counter" . nfIO $ incrsIO $ counters   0.1 0.1
--      , bench "eps = 0.1 d = 0.1 counter'" . nfIO $ incrsIO $ counters' 0.1 0.1

--      , bench "eps = 0.05 d = 0.05 counter" . nfIO $ incrsIO $ counters   0.05 0.05
--      , bench "eps = 0.05 d = 0.05 counter'" . nfIO $ incrsIO $ counters' 0.05 0.05
--  ]

--incrsIO :: MonadRandom m => [[Counter]] -> m [[Counter]]
--incrsIO  xxs = runRVar (incrs  xxs) StdRandom

--counters :: Eps -> Delta  -> [[Counter]]
--counters e d = replicate m $ replicate t 0
--  where (m,t) = (m' e d, t' e d)

--counters' :: Eps -> Delta -> [[Counter]]
--counters' e d = replicate t $ replicate m 0
--  where (m,t) = (m' e d, t' e d)

--t' e d = round $ 1/(e^2*d)
--m' e d = round . log $ 1/d

---- * given a list of list of counters toss a coin for each counter and incr
--incrs :: [[Counter]] -> RVar [[Counter]]
--incrs = sequence . fmap (sequence . fmap incr)

---- * Increment a counter `x` with probability 1/2^x
--incr :: Counter -> RVar Counter
--incr x = do
--  h <- toss . coin $ 0.5^(round x)
--  return $ if isHead h then (seq () succ x) else seq () x

{-----------------------------------------------------------------------------
    Benchmark time compleity of mean/median
       see : http://www.serpentine.com/criterion/tutorial.html
------------------------------------------------------------------------------}

-- * observation: 5% difference in runtime using mean   vs not
-- *              5% difference in runtime using median vs not
--main :: IO ()
--main = defaultMain . return $ bgroup "morris mean vs no median" [
--      bench "morris median    100 items eps = 0.1   delta = 0.1 " . nfIO $ morris  0.1 0.1   [1..100]
--    , bench "morris no median 100 items eps = 0.1   delta = 0.1 " . nfIO $ morris' 0.1 0.1   [1..100]
--    , bench "morris    median 100 items eps = 0.05  delta = 0.05" . nfIO $ morris  0.05 0.05 [1..100]
--    , bench "morris no median 100 items eps = 0.05  delta = 0.05" . nfIO $ morris' 0.05 0.05 [1..100]
--  ]

{-----------------------------------------------------------------------------
    II. Approximate Mean
------------------------------------------------------------------------------}

-- * Run Morris alpha on stream inputs `xs`
--morrisA :: [a] -> IO Counter
--morrisA xs = flip runRVar StdRandom $ Cl.sourceList xs $$ alpha

---- * Run Morris beta on stream inputs `xs` for `t` independent trials and average
--morrisB :: Int -> [a] -> IO Counter
--morrisB t =  fmap rmean . replicateM t . morrisA

---- * final morris algorithm
---- * Run on stream inputs `xs` for t independent trials for `t = 1/eps`,
---- * and `m` times in parralell, for `m = 1/(e^2 * d)`
---- * and take the median
---- * TODO: make this actually parralell
--morris :: Eps -> Delta -> [a] -> IO Counter
--morris e d = fmap rmedian . replicateM m . morrisB t
--  where (t,m) = (round $ 1/(e^2*d), round $ 1/d)

------------------------------------------------------------------------------
--    III. Utils
-------------------------------------------------------------------------------

---- * Utils * --

---- * A step in morris Algorithm alpha
--alpha :: Sink a RVar Counter
--alpha = (\x -> 2^(round x) - 1) <$> Cl.foldM (\x _ -> incr x) 0

---- * Increment a counter `x` with probability 1/2^x
--incr :: Counter -> RVar Counter
--incr x = do
--  h <- toss . coin $ 0.5^(round x)
--  return $ if isHead h then (seq () succ x) else seq () x

--rmean, rmedian :: (Floating a, Ord a, RealFrac a) => [a] -> Float
--rmean   = fromIntegral . round . mean
--rmedian = fromIntegral . round . median

{-----------------------------------------------------------------------------
    Depricated
------------------------------------------------------------------------------}

-- * Non fmap f version to test list traversal
-- * marginal difference due to fmap g where g = mean or median
--morrisB' :: Int -> [a] -> IO [Counter]
--morrisB' t = replicateM t . morrisA

--morris' :: Eps -> Delta -> [a] -> IO [Counter]
--morris' e d = replicateM m . morrisB t
  --where (t,m) = (round $ 1/(e^2*d), round $ 1/d)

{-----------------------------------------------------------------------------
  Benchmark  median of means
  naive "imperitive" solution worse in large n
------------------------------------------------------------------------------}

--main :: IO ()
--main = defaultMain . return $ bgroup "count vs count'" [
--      --  bench "median  1000"  $ whnf (medianOfMeans  10) [1..1000]
--      --, bench "median' 1000"  $ whnf (medianOfMeans' 10) [1..1000]

--      --, bench "median 10000"  $ whnf (medianOfMeans  10) [1..10000]
--      --, bench "median' 10000"  $ whnf (medianOfMeans' 10) [1..10000]

--      --, bench "median 100000"  $ whnf (medianOfMeans  10) [1..100000]
--      --, bench "median' 100000"  $ whnf (medianOfMeans' 10) [1..100000]

--        bench "median  0.1"  $ whnf (medianOfMeans  1000000) [1..5000000]
--      , bench "median' 0.1"  $ whnf (medianOfMeans' 1000000) [1..5000000]

--      , bench "median  0.5"  $ whnf (medianOfMeans  8000) [1..24000]
--      , bench "median' 0.5"  $ whnf (medianOfMeans' 8000) [1..24000]

--  ]

---- * this is worse!!
--medianOfMeans' :: Counter -> [Counter] -> Counter
--medianOfMeans' t xs = let (ms,_,_) = foldr (tomean t) ([],0,1) xs in median ms

--tomean :: Counter -> Counter -> ([Counter],Counter,Counter) -> ([Counter],Counter,Counter)
--tomean t x (ms,m,c) | c < t      = (ms,m+x,c+1)
--                    | otherwise  = ((fromIntegral . round $ (m+x)/t):ms,0,1)

--medianOfMeans :: Int -> [Counter] -> Counter
--medianOfMeans t = median . fmap mean' . (chunksOf t)
-- fmap median' .

--ms = medianOfMeans 10 [1..100]
--(ms',_,_) = medianOfMeans' 10 [1..100]


text:
{-# LANGUAGE StandaloneDeriving #-}
module Examples.Quantum where

import LAoP.Utils
import LAoP.Matrix.Type
import Data.Complex
import Prelude hiding (id, (.))

deriving instance Ord a => Ord (Complex a)

xor :: (Bool, Bool) -> Bool
xor (False, b) = b
xor (True, b) = not b

ker :: Num e => Matrix e a b -> Matrix e a a
ker m = tr m . m

cnot :: Matrix (Complex Double) (Bool, Bool) (Bool, Bool)
cnot = kr fstM (fromF xor)

ccnot :: (Num e, Ord e) => Matrix e ((Bool, Bool), Bool) ((Bool, Bool), Bool)
ccnot = kr fstM (fromF f)
  where
    f = xor . tp (uncurry (&&)) id
    tp f g (a,b) = (f a, g b)

had :: Matrix (Complex Double) Bool Bool
had = (1/sqrt 2) .| fromLists [[1, 1], [1, -1]]

bell :: Matrix (Complex Double) (Bool, Bool) (Bool, Bool)
bell = cnot . (had >< iden)


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE InstanceSigs #-}
{-# LANGUAGE Rank2Types #-}
{-# LANGUAGE FlexibleInstances    #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE IncoherentInstances  #-}
{-# LANGUAGE StandaloneDeriving  #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ForeignFunctionInterface #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

module Data.Vector.Image.Color.RGB (
  RGB(..)
, ColorSpace(..)
, RGB8
) where

import Prelude hiding(map)
import qualified Prelude as P
import qualified Data.List as L
import Data.Vector.Storable hiding(forM_,map,replicate,fromList,toList,(!?),(!))
import Data.Maybe
import Numeric.LinearAlgebra.Data hiding (fromList,toList,(!))
import Numeric.LinearAlgebra.HMatrix hiding (fromList,toList,(!))
import qualified Data.Vector.Storable as V
import Foreign
import GHC.Real
import Foreign.C.String
import Foreign.C.Types
import qualified Control.Monad as C

data RGB a = RGB {
   rgb_r :: a
 , rgb_g :: a
 , rgb_b :: a
} deriving (Show,Eq,Ord)

type RGB8 = RGB Word8

class ColorSpace color b where
  toRGB :: color b -> RGB b
  fromRGB :: RGB b -> color b
  toList :: color b -> [b]
  fromList :: [b] -> color b
  dim :: color b -> Int
  dummyElem :: color b -> b
  dummyElem _ = error "do not eval this"

instance Functor RGB where
  fmap func (RGB r g b) = RGB (func r) (func g) (func b)

op :: (ColorSpace c a,Num a) => (a -> a -> a) -> c a -> c a -> c a
op func rgb rgb2 =
  let rgb' = toList rgb
      rgb2' = toList rgb2
      in fromList $ P.map (uncurry func) $ zip rgb' rgb2'

fmap' :: (ColorSpace c a, ColorSpace c b)
     => (a -> b)
     -> c a
     -> c b
fmap' func val = fromList $ P.map func $ toList val

instance (ColorSpace c a,Num a) => Num (c a) where
  (+) = op (+)
  (-) = op (-)
  (*) = op (*)
  negate = fmap' negate
  abs = fmap' abs
  signum = fmap' signum
  fromInteger v = fromList $ replicate (dim dummyColor) v'
    where
      v' = fromIntegral v
      dummyColor = error "do not eval this" :: (c a)

instance Enum (RGB Int) where
  succ rgb =  rgb + RGB 1 1 1
  pred rgb =  rgb - RGB 1 1 1
  toEnum i = fromIntegral i
  fromEnum (RGB r g b) = r
  enumFrom rgb = rgb : (enumFrom $ succ rgb)
  enumFromThen rgb rgb1 = rgb : (enumFromThen rgb1 (2 * rgb1 - rgb))

instance Real (RGB Int) where
  toRational (RGB r g b) = (fromIntegral $ (r + 2*g + b) `div` 4) % 1

instance Integral (RGB Int) where
  quot = op quot
  rem = op rem
  div = op div
  mod = op mod
  quotRem a b = (quot a b,rem a b)
  divMod a b = (div a b,mod a b)
  toInteger (RGB r g b) = fromIntegral $ (r + 2*g + b) `div` 4

instance (Storable a) => Element a
instance (ColorSpace c a,Storable a) => Element (c a)

instance (ColorSpace c a,Storable a) => Storable (c a) where
  sizeOf    d = dim d * sizeOf (dummyElem d)
  alignment d = dim d * sizeOf (dummyElem d)
  peek ptr = do
    v <- C.forM [0..((dim d)-1)] $ \idx -> do
      peekByteOff ptr (sizeOf (dummyElem d) * idx)
    return $ fromList v
    where
      d = error "do not eval this" :: c a
  poke ptr val = do
    C.forM_ (zip [0..] (toList val)) $ \(idx,v) -> do
      pokeByteOff ptr (sizeOf (dummyElem d) * idx) v
    where
      d = error "do not eval this" :: c a

instance ColorSpace RGB a where
  toRGB = id
  fromRGB = id
  toList (RGB r g b) = [r,g,b]
  fromList [r,g,b] = (RGB r g b)
  fromList val = error $ "List length must be 3 but " L.++  show (L.length val)
  dim _ = 3

max_stars_repo_name: junjihashimoto/img-vector

max_issues_repo_name: junjihashimoto/img-vector

max_forks_repo_name: junjihashimoto/img-vector


text:
{-# LANGUAGE TypeApplications #-}
{-# OPTIONS_GHC -fno-warn-type-defaults #-}

module ScalarWave where

import Control.Exception (assert)
import Control.Monad.Loops
import qualified Data.Vector.Unboxed as U
import Data.VectorSpace
import Dual
import Equations
import Foreign
import qualified Numeric.LinearAlgebra as L
import SSSFN
import System.CPUTime

default (Int)

--------------------------------------------------------------------------------

bndw :: Floating a => Coord a -> a
bndw (Coord u v) = swamp * cos (pi / 4 * swk * t) * cos (pi / 4 * swk * r)
  where
    t = v + u
    r = v - u
    swk = 0.1
    swamp = 0.1

guess1 :: Floating a => Coord a -> State a
guess1 c = State 1 1 0

guess :: (L.Container L.Vector a, Floating a) => Grid (State a)
guess = gsample (\(u, v) -> guess1 (Coord u v))

--------------------------------------------------------------------------------

eqnInt1 :: Num a => Eqn -> Coord a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> a
eqnInt1 eqn c s r2s r1sr r1su r1sv su sv suv suu svv =
  case eqn of
    EqnG01 -> eqnG01 s r2s r1sr su sv suv
    EqnG22 -> eqnG22 s r2s r1sr su sv suv
    EqnSW2 -> eqnSW2 s r2s r1sr su sv suv
    EqnG00 -> eqnG00 c s r2s r1su su suu
    EqnG11 -> eqnG11 c s r2s r1sv sv svv

eqnOrig1 :: Floating a => Eqn -> Coord a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> a
eqnOrig1 eqn c s r2s r1sr r1su r1sv su sv suv suu svv =
  case eqn of
    EqnG01 -> q s - 1
    EqnG22 -> g s - 1
    EqnSW2 -> w s - bndw c
    EqnG00 -> 0
    EqnG11 -> 0

-- i == 0, u == -1
eqnBndU1 :: (Eq a, Floating a) => Eqn -> Coord a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> a
eqnBndU1 eqn c s r2s r1sr r1su r1sv su sv suv suu svv =
  assert (u c == -1) $
    case eqn of
      EqnG01 -> eqnG11 c s r2s r1sv sv svv
      EqnG22 -> g s - 1
      EqnSW2 -> w s - bndw c
      EqnG00 -> 0
      EqnG11 -> 0

-- j == 0, v == -1
eqnBndV1 :: (Eq a, Floating a) => Eqn -> Coord a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> a
eqnBndV1 eqn c s r2s r1sr r1su r1sv su sv suv suu svv =
  assert (v c == -1) $
    case eqn of
      EqnG01 -> eqnG00 c s r2s r1su su suu
      EqnG22 -> g s - 1
      EqnSW2 -> w s - bndw c
      EqnG00 -> 0
      EqnG11 -> 0

eqnEndU1 :: Floating a => Eqn -> Coord a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> a
eqnEndU1 eqn c s r2s r1sr r1su r1sv su sv suv suu svv =
  case eqn of
    EqnG01 -> q s - 1
    EqnG22 -> g s - 1
    EqnSW2 -> w s - bndw c
    EqnG00 -> 0
    EqnG11 -> 0

eqnEndV1 :: Floating a => Eqn -> Coord a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> a
eqnEndV1 eqn c s r2s r1sr r1su r1sv su sv suv suu svv =
  case eqn of
    EqnG01 -> q s - 1
    EqnG22 -> g s - 1
    EqnSW2 -> w s - bndw c
    EqnG00 -> 0
    EqnG11 -> 0

equation1 :: (Eq a, Floating a) => Coord Int -> Eqn -> Coord a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> State a -> a
equation1 (Coord i j) =
  if
    | i == 0 && j == 0 -> eqnOrig1
    | i == 0 && j == np - 1 -> eqnEndU1
    | i == np - 1 && j == 0 -> eqnEndV1
    | i == 0 -> eqnBndU1
    | j == 0 -> eqnBndV1
    | True -> eqnInt1

--------------------------------------------------------------------------------

g2s :: Storable a => Grid a -> State (Grid a)
g2s (Grid xs) =
  let [g, q, w] = L.takesV (replicate 3 (np * np)) xs
   in State (Grid g) (Grid q) (Grid w)

s2g :: Storable a => State (Grid a) -> Grid a
s2g (State (Grid g) (Grid q) (Grid w)) = Grid (L.vjoin [g, q, w])

eqn1 ::
  forall a.
  (L.Field a, RealFloat a, U.Unbox a) =>
  Eqn ->
  Coord (Grid a) ->
  State (Grid a) ->
  Grid a
eqn1 eqn c s =
  Grid $
    L.fromList
      [ eqnRow i j
        | i <- [0 .. np - 1],
          j <- [0 .. np - 1]
      ]
  where
    r2s = (r2project #>) <$> s
    r1sr = (r1derivR #>) <$> s
    r1su = (r1derivU #>) <$> s
    r1sv = (r1derivV #>) <$> s
    su = (derivU #>) <$> s
    sv = (derivV #>) <$> s
    suv = (derivUV #>) <$> s
    suu = (derivUU #>) <$> s
    svv = (derivVV #>) <$> s
    eqnRow :: Int -> Int -> a
    eqnRow i j =
      let at :: Grid a -> a
          at x = x ! (i, j)
       in equation1
            (Coord i j)
            eqn
            (at <$> c)
            (at <$> s)
            (at <$> r2s)
            (at <$> r1sr)
            (at <$> r1su)
            (at <$> r1sv)
            (at <$> su)
            (at <$> sv)
            (at <$> suv)
            (at <$> suu)
            (at <$> svv)

eqns ::
  forall a.
  (L.Field a, RealFloat a, U.Unbox a) =>
  Coord (Grid a) ->
  State (Grid a) ->
  Grid a
eqns c s =
  Grid $
    L.fromList
      [ eqnRow (i, j, eqn)
        | eqn <- [EqnG01, EqnG22, EqnSW2],
          i <- [0 .. np - 1],
          j <- [0 .. np - 1]
      ]
  where
    r2s = (r2project #>) <$> s
    r1sr = (r1derivR #>) <$> s
    r1su = (r1derivU #>) <$> s
    r1sv = (r1derivV #>) <$> s
    su = (derivU #>) <$> s
    sv = (derivV #>) <$> s
    suv = (derivUV #>) <$> s
    suu = (derivUU #>) <$> s
    svv = (derivVV #>) <$> s
    eqnRow :: (Int, Int, Eqn) -> a
    eqnRow (i, j, eqn) =
      let at :: Grid a -> a
          at x = x ! (i, j)
       in equation1
            (Coord i j)
            eqn
            (at <$> c)
            (at <$> s)
            (at <$> r2s)
            (at <$> r1sr)
            (at <$> r1su)
            (at <$> r1sv)
            (at <$> su)
            (at <$> sv)
            (at <$> suv)
            (at <$> suu)
            (at <$> svv)

eqns5 ::
  forall a.
  (L.Field a, RealFloat a, U.Unbox a) =>
  Coord (Grid a) ->
  State (Grid a) ->
  Grid a
eqns5 c s =
  Grid $
    L.fromList
      [ eqnRow (i, j, eqn)
        | eqn <- [EqnG01, EqnG22, EqnSW2, EqnG00, EqnG11],
          i <- [0 .. np - 1],
          j <- [0 .. np - 1]
      ]
  where
    r2s = (r2project #>) <$> s
    r1sr = (r1derivR #>) <$> s
    r1su = (r1derivU #>) <$> s
    r1sv = (r1derivV #>) <$> s
    su = (derivU #>) <$> s
    sv = (derivV #>) <$> s
    suv = (derivUV #>) <$> s
    suu = (derivUU #>) <$> s
    svv = (derivVV #>) <$> s
    eqnRow :: (Int, Int, Eqn) -> a
    eqnRow (i, j, eqn) =
      let at :: Grid a -> a
          at x = x ! (i, j)
       in equation1
            (Coord i j)
            eqn
            (at <$> c)
            (at <$> s)
            (at <$> r2s)
            (at <$> r1sr)
            (at <$> r1su)
            (at <$> r1sv)
            (at <$> su)
            (at <$> sv)
            (at <$> suv)
            (at <$> suu)
            (at <$> svv)

op ::
  forall a.
  (L.Container L.Vector a, L.Field a, RealFloat a, U.Unbox a) =>
  Coord (Grid a) ->
  State (Grid a) ->
  Op a
op (Coord u v) s =
  Op $
    L.fromRows
      [ jacRow (i, j, eqn)
        | eqn <- [EqnG01, EqnG22, EqnSW2],
          i <- [0 .. np - 1],
          j <- [0 .. np - 1]
      ]
  where
    r2s = (r2project #>) <$> s
    r1sr = (r1derivR #>) <$> s
    r1su = (r1derivU #>) <$> s
    r1sv = (r1derivV #>) <$> s
    su = (derivU #>) <$> s
    sv = (derivV #>) <$> s
    suv = (derivUV #>) <$> s
    suu = (derivUU #>) <$> s
    svv = (derivVV #>) <$> s
    jacRow :: (Int, Int, Eqn) -> L.Vector a
    jacRow (i, j, eqn) =
      L.vjoin $
        [ let con :: Grid a -> Dual a
              der :: Grid a -> Dual a
              con x = Dual (x ! (i, j)) 0
              der x = Dual (x ! (i, j)) 1
              scon :: State (Grid a) -> State (Dual a)
              sder :: State (Grid a) -> State (Dual a)
              scon (State g q w) = State (con g) (con q) (con w)
              sder (State g q w) = case var of
                VarG -> State (der g) (con q) (con w)
                VarQ -> State (con g) (der q) (con w)
                VarW -> State (con g) (con q) (der w)
           in ( ( dual $
                    equation1
                      (Coord i j)
                      eqn
                      (Coord (con u) (con v))
                      (sder s)
                      (scon r2s)
                      (scon r1sr)
                      (scon r1su)
                      (scon r1sv)
                      (scon su)
                      (scon sv)
                      (scon suv)
                      (scon suu)
                      (scon svv)
                )
                  `L.scale` (getOp delta L.! (i * np + j))
              )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (sder r2s)
                                (scon r1sr)
                                (scon r1su)
                                (scon r1sv)
                                (scon su)
                                (scon sv)
                                (scon suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp r2project L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon r2s)
                                (sder r1sr)
                                (scon r1su)
                                (scon r1sv)
                                (scon su)
                                (scon sv)
                                (scon suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp r1derivR L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon r2s)
                                (scon r1sr)
                                (sder r1su)
                                (scon r1sv)
                                (scon su)
                                (scon sv)
                                (scon suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp r1derivU L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon r2s)
                                (scon r1sr)
                                (scon r1su)
                                (sder r1sv)
                                (scon su)
                                (scon sv)
                                (scon suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp r1derivV L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon r2s)
                                (scon r1sr)
                                (scon r1su)
                                (scon r1sv)
                                (sder su)
                                (scon sv)
                                (scon suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp derivU L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon r2s)
                                (scon r1sr)
                                (scon r1su)
                                (scon r1sv)
                                (scon su)
                                (sder sv)
                                (scon suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp derivV L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon r2s)
                                (scon r1sr)
                                (scon r1su)
                                (scon r1sv)
                                (scon su)
                                (scon sv)
                                (sder suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp derivUV L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon r2s)
                                (scon r1sr)
                                (scon r1su)
                                (scon r1sv)
                                (scon su)
                                (scon sv)
                                (scon suv)
                                (sder suu)
                                (scon svv)
                          )
                            `L.scale` (getOp derivUU L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon r2s)
                                (scon r1sr)
                                (scon r1su)
                                (scon r1sv)
                                (scon su)
                                (scon sv)
                                (scon suv)
                                (scon suu)
                                (sder svv)
                          )
                            `L.scale` (getOp derivVV L.! (i * np + j))
                        )
          | var <- [VarG, VarQ, VarW]
        ]

op5 ::
  forall a.
  (L.Container L.Vector a, L.Field a, RealFloat a, U.Unbox a) =>
  Coord (Grid a) ->
  State (Grid a) ->
  Op a
op5 (Coord u v) s =
  Op $
    L.fromRows
      [ jacRow (i, j, eqn)
        | eqn <- [EqnG01, EqnG22, EqnSW2, EqnG00, EqnG11],
          i <- [0 .. np - 1],
          j <- [0 .. np - 1]
      ]
  where
    r2s = (r2project #>) <$> s
    r1sr = (r1derivR #>) <$> s
    r1su = (r1derivU #>) <$> s
    r1sv = (r1derivV #>) <$> s
    su = (derivU #>) <$> s
    sv = (derivV #>) <$> s
    suv = (derivUV #>) <$> s
    suu = (derivUU #>) <$> s
    svv = (derivVV #>) <$> s
    jacRow :: (Int, Int, Eqn) -> L.Vector a
    jacRow (i, j, eqn) =
      L.vjoin $
        [ let con :: Grid a -> Dual a
              der :: Grid a -> Dual a
              con x = Dual (x ! (i, j)) 0
              der x = Dual (x ! (i, j)) 1
              scon :: State (Grid a) -> State (Dual a)
              sder :: State (Grid a) -> State (Dual a)
              scon (State g q w) = State (con g) (con q) (con w)
              sder (State g q w) = case var of
                VarG -> State (der g) (con q) (con w)
                VarQ -> State (con g) (der q) (con w)
                VarW -> State (con g) (con q) (der w)
           in ( ( dual $
                    equation1
                      (Coord i j)
                      eqn
                      (Coord (con u) (con v))
                      (sder s)
                      (scon r2s)
                      (scon r1sr)
                      (scon r1su)
                      (scon r1sv)
                      (scon su)
                      (scon sv)
                      (scon suv)
                      (scon suu)
                      (scon svv)
                )
                  `L.scale` (getOp delta L.! (i * np + j))
              )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (sder r2s)
                                (scon r1sr)
                                (scon r1su)
                                (scon r1sv)
                                (scon su)
                                (scon sv)
                                (scon suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp r2project L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon r2s)
                                (sder r1sr)
                                (scon r1su)
                                (scon r1sv)
                                (scon su)
                                (scon sv)
                                (scon suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp r1derivR L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon r2s)
                                (scon r1sr)
                                (sder r1su)
                                (scon r1sv)
                                (scon su)
                                (scon sv)
                                (scon suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp r1derivU L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon r2s)
                                (scon r1sr)
                                (scon r1su)
                                (sder r1sv)
                                (scon su)
                                (scon sv)
                                (scon suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp r1derivV L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon r2s)
                                (scon r1sr)
                                (scon r1su)
                                (scon r1sv)
                                (sder su)
                                (scon sv)
                                (scon suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp derivU L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon r2s)
                                (scon r1sr)
                                (scon r1su)
                                (scon r1sv)
                                (scon su)
                                (sder sv)
                                (scon suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp derivV L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon r2s)
                                (scon r1sr)
                                (scon r1su)
                                (scon r1sv)
                                (scon su)
                                (scon sv)
                                (sder suv)
                                (scon suu)
                                (scon svv)
                          )
                            `L.scale` (getOp derivUV L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon r2s)
                                (scon r1sr)
                                (scon r1su)
                                (scon r1sv)
                                (scon su)
                                (scon sv)
                                (scon suv)
                                (sder suu)
                                (scon svv)
                          )
                            `L.scale` (getOp derivUU L.! (i * np + j))
                        )
                `L.add` ( ( dual $
                              equation1
                                (Coord i j)
                                eqn
                                (Coord (con u) (con v))
                                (scon s)
                                (scon r2s)
                                (scon r1sr)
                                (scon r1su)
                                (scon r1sv)
                                (scon su)
                                (scon sv)
                                (scon suv)
                                (scon suu)
                                (sder svv)
                          )
                            `L.scale` (getOp derivVV L.! (i * np + j))
                        )
          | var <- [VarG, VarQ, VarW]
        ]

--------------------------------------------------------------------------------

getTime :: IO Double
getTime = do
  t <- getCPUTime
  return $ fromInteger t / 1.0e+12

main0 :: IO ()
main0 =
  do
    putStrLn $ show $ gmaxabs $ coordU @Double
    putStrLn $ show $ gmaxabs $ coordV @Double
    putStrLn $ show $ gmaxabs $ coordR @Double
    putStrLn $ show $ gmaxabs $ coordR1 @Double
    putStrLn $ show $ omaxabs $ zproject @Double
    putStrLn $ show $ omaxabs $ r2project @Double
    putStrLn $ show $ omaxabs $ r1derivR @Double
    putStrLn $ show $ omaxabs $ r1derivU @Double
    putStrLn $ show $ omaxabs $ r1derivV @Double
    putStrLn $ show $ omaxabs $ derivU @Double
    putStrLn $ show $ omaxabs $ derivV @Double
    putStrLn $ show $ omaxabs $ derivUV @Double
    putStrLn $ show $ omaxabs $ derivUU @Double
    putStrLn $ show $ omaxabs $ derivVV @Double

main :: IO ()
main =
  do
    putStrLn "Spherically Symmetric Scalar Field in double-Null coordinates"
    let coords = Coord coordU (coordV @Double)
    let ini = State (gmap g guess) (gmap q guess) (gmap w guess)
    vars0 <- return (s2g ini)
    let done (n, vars) =
          let res = eqns coords (g2s vars)
              maxabsres = gmaxabs res
           in isNaN maxabsres
                || isInfinite maxabsres
                || maxabsres < 1.0e-12
                || n >= 5
    let iter (n, vars) = do
          let res = eqns coords (g2s vars)
          let res5 = eqns5 coords (g2s vars)
          t <- getTime
          putStrLn $
            ("iter " ++ show n)
              ++ ("   time " ++ show (round t) ++ " sec")
              ++ ("   |res| " ++ show (gmaxabs res))
              ++ ("   |res5| " ++ show (gmaxabs res5))
          let jac = op coords (g2s vars)
          let dvars = solve jac res
          return (n + 1, vars ^-^ dvars)
    let done5 (n, vars) =
          let res5 = eqns5 coords (g2s vars)
              maxabsres5 = gmaxabs res5
           in isNaN maxabsres5
                || isInfinite maxabsres5
                || maxabsres5 < 1.0e-12
                || n >= 5
    let iter5 (n, vars) = do
          let res = eqns coords (g2s vars)
          let res5 = eqns5 coords (g2s vars)
          t <- getTime
          putStrLn $
            ("iter " ++ show n)
              ++ ("   time " ++ show (round t) ++ " sec")
              ++ ("   |res| " ++ show (gmaxabs res))
              ++ ("   |res5| " ++ show (gmaxabs res5))
          let jac5 = op5 coords (g2s vars)
          let dvars = solveLS jac5 res5
          return (n + 1, vars ^-^ dvars)
    -- (n, vars) <- iterateUntilM done iter (0, vars0)
    (n, vars) <- iterateUntilM done5 iter5 (0, vars0)
    let res = eqns coords (g2s vars)
    let res5 = eqns5 coords (g2s vars)
    t <- getTime
    putStrLn $
      ("iter " ++ show n)
        ++ ("   time " ++ show (round t) ++ " sec")
        ++ ("   |res| " ++ show (gmaxabs res))
        ++ ("   |res5| " ++ show (gmaxabs res5))
    putStrLn $ "var.g " ++ show (gmap approx $ g $ g2s vars)
    putStrLn $ "var.q " ++ show (gmap approx $ q $ g2s vars)
    putStrLn $ "var.w " ++ show (gmap approx $ w $ g2s vars)
    putStrLn $ "res.G01 " ++ show (gmap approx $ eqn1 EqnG01 coords (g2s vars))
    putStrLn $ "res.G22 " ++ show (gmap approx $ eqn1 EqnG22 coords (g2s vars))
    putStrLn $ "res.SW2 " ++ show (gmap approx $ eqn1 EqnSW2 coords (g2s vars))
    putStrLn $ "res.G00 " ++ show (gmap approx $ eqn1 EqnG00 coords (g2s vars))
    putStrLn $ "res.G11 " ++ show (gmap approx $ eqn1 EqnG11 coords (g2s vars))
    putStrLn $ "S var.g " ++ show (gmap approx $ inv vandermonde #> (g $ g2s vars))
    putStrLn $ "S var.q " ++ show (gmap approx $ inv vandermonde #> (q $ g2s vars))
    putStrLn $ "S var.w " ++ show (gmap approx $ inv vandermonde #> (w $ g2s vars))
    putStrLn $ "S res.G01 " ++ show (gmap approx $ inv vandermonde #> (eqn1 EqnG01 coords (g2s vars)))
    putStrLn $ "S res.G22 " ++ show (gmap approx $ inv vandermonde #> (eqn1 EqnG22 coords (g2s vars)))
    putStrLn $ "S res.SW2 " ++ show (gmap approx $ inv vandermonde #> (eqn1 EqnSW2 coords (g2s vars)))
    putStrLn $ "S res.G00 " ++ show (gmap approx $ inv vandermonde #> (eqn1 EqnG00 coords (g2s vars)))
    putStrLn $ "S res.G11 " ++ show (gmap approx $ inv vandermonde #> (eqn1 EqnG11 coords (g2s vars)))
    putStrLn "Done."
  where
    approx :: RealFrac a => a -> a
    approx x = fromInteger (round (1.0e+10 * x)) / 1.0e+10

max_stars_repo_name: eschnett/sssfn

max_issues_repo_name: eschnett/sssfn

max_forks_repo_name: eschnett/sssfn


text:
module Lib
    ( euclidianDistance
    , Point
    , Vector
    , makePoint
    , eye1
    , getEyeCenter
    , getLineFromPoints
    , projection
    ) where

import qualified Numeric.LinearAlgebra.HMatrix as HM

type Point = HM.Vector Double
type Vector = HM.Vector Double
data Line = Line { lineCoefficientMatrix :: HM.Matrix Double, lineTerms :: HM.Vector Double } deriving Show

makePoint :: Double -> Double -> Double -> Double -> Point
makePoint x y z t = HM.vector [x, y, z, t]

makeVector :: Double -> Double -> Double -> Double -> Vector
makeVector = makePoint

data EyeConfig = EyeConfig { retinaRadius :: Double, convergenceDistance :: Double }

data EyePosition = EyePosition { convergencePoint :: Point, direction :: Vector }

data Eye = Eye { config :: EyeConfig, position :: EyePosition }

eye1 :: Eye
eye1 = Eye {
    config = EyeConfig { retinaRadius = 10, convergenceDistance = 30 },
    position = EyePosition { convergencePoint = makePoint 1 1 1 1, direction = makeVector (-1) (-1) (-1) (-1) }
    }

euclidianDistance :: Point -> Point -> Double
euclidianDistance a b = sqrt . sum . map (**2) $ zipWith (-) (HM.toList a) (HM.toList b)

euclidianLength :: Vector -> Double
euclidianLength = euclidianDistance (makePoint 0 0 0 0)

getEyeCenter :: Eye -> Point
getEyeCenter Eye { config = eyeConfig, position = eyePosition } =
  let
    ratio = (/ (euclidianLength . direction $ eyePosition)) (convergenceDistance eyeConfig)
  in
    (+ convergencePoint eyePosition) . (* direction eyePosition) . HM.scalar $ ratio

getLineFromPoints :: Point -> Point -> Line
getLineFromPoints a b =
  let
    l = a - b
    a' = HM.toList a
    l' = HM.toList l
    coefficientXy = [l'!!1, -(head l'), 0, 0]
    coefficientYz = [0, l'!!2, -(l'!!1), 0]
    coefficientZt = [0, 0, l'!!3, -(l'!!2)]
    termXy = HM.vector coefficientXy HM.<.> a
    termYz = HM.vector coefficientYz HM.<.> a
    termZt = HM.vector coefficientZt HM.<.> a
  in
    Line {
        lineCoefficientMatrix = HM.matrix 4 (coefficientXy ++ coefficientYz ++ coefficientZt),
        lineTerms = HM.vector [termXy, termYz, termZt]
    }

projection :: Eye -> Point -> Maybe Point
projection eye point =
  let
    Eye { config = eyeConfig, position = eyePosition } = eye
    eyeCenter = getEyeCenter eye
    Line { lineCoefficientMatrix = coefficients, lineTerms = terms } =
        getLineFromPoints point (convergencePoint eyePosition)
    eyeDirection = direction eyePosition
    intersectionCoefficients = HM.fromRows . (++ [eyeDirection]) . HM.toRows $ coefficients
    intersectionTerms = HM.matrix 1 . (++ [eyeDirection HM.<.> eyeCenter]) . HM.toList $ terms
  in
    HM.flatten <$> HM.linearSolve intersectionCoefficients intersectionTerms


text:
module Lec7 where

-- Функторы
--
-- Функция map была бы логична не только для списков. Как бы она работала для Maybe?
-- f x = x + 1
-- x1 = Just 42
-- x2 = Nothing
--
-- -хочу применить f к x1, x2, точнее, к значению внутри
-- maybeMap fun (Just n) = Just (fun n)
-- maybeMap _ Nothing = Nothing
-- print $ maybeMap f x1  -> Just 43
-- print $ maybeMap f x2  -> Nothing
--
-- Допустим, я хочу уметь так же применять функции ко вторым элементам пары:
-- f x = x + 1
-- x1 = ("hello", 42)
-- x2 = ([1, 2, 3], 42.0)
-- хочу применить f ко второму элементу
-- pairMap fun (x, y) = (x, fun y)
-- print $ pairMap f x1  -> ("hello", 43)
-- print $ pairMap f x2  -> ([1, 2, 3], 43.0)
--
-- Есть и другие аналогичные ситуации. Когда значение (или значения) хранятся в какой-то структуре с дополнительными
-- данными.
-- Список - несколько значений и информация об их порядке
-- Maybe - одно значение или ноль
-- Пара - значение вместе с дополнительной информацией
--
-- Функции тоже можно понимать как такие структуры:
-- fun :: a -> b
-- получается, функция хранит значение типа B, доступ к которым осуществляется по аргументу типа a:
-- Например,
-- f 1 = "один"
-- f 2 = "два"
-- функция "хранит" внутри себя значения "один", "два", получить их можно, подставив соответствующий аргумент в функцию.
-- f x = x + 1
-- x1 = \x -> 2 * x
-- x2 = \x -> 3 * x
-- funMap f x1 - должна получиться функция, которая сопадает с x1, но после применения x1 дополнительно вызывает f
-- funMap fun x = \t -> fun (x t)
--
-- xx1 = funMap f x1
-- xx2 = funMap f x2
-- print $ xx1 42  -> 85
-- print $ xx2 42  -> 127
--
-- На самом деле все эти map функции уже есть. Они называются fmap или <$>, и они делают ровно то же, что делаем мы:
-- print $ fmap (+1) (Just 42)  -> Just 43
-- print $ (+1) <$> Just 42  -> Just 43
-- print $ (+1) <$> [10, 20, 30]  -> [11, 21, 31]
-- print $ (+1) <$> ("abc", 123)  -? ("abc", 124)
--
-- Это работает, потому что в Haskell есть класс типов Functor - функтор. Он определен так:
-- class Functor f where
--   fmap :: (a -> b) -> f a -> f b
--
-- У него есть реализации, например, примерно так выглядит реализация для Maybe.
-- instance Functor Maybe where
--   fmap f (Just x) = Just (f x)
--   fmap _ Nothing = Nothing
--
-- Если вы определяете функтор, вы должны гарантировать выполнение следующего правила. Функция id - это функция, которая
-- ничего не делает: id x = x, её тип a -> a
-- **Правило**:
-- 1) fmap id x == x (или эквивалентно) fmap id == id
-- 2) fmap (g . f) x == fmap g . fmap f
--
-- ## Аппликативные функторы
-- Расширение обычных функторов дополнительными операциями. Представьте опять же Maybe.
-- Было в функторах:
-- 1) Обычная функция
-- 2) Данные внутри структуры
-- 3) Применяем функцию к каждому данному внутри структуры
-- Теперь в аппликативных функторах можно дополнительно:
-- 1) Обычные функция внутри структуры
-- 2) Данные внутри структуры
-- 3) Применяем функции из структуры к данным. Как именно применять - очень сильно зависит от ситуации.
--
-- Пример с Maybe:
-- f = Just (+1) -- увеличение на 1, но внутри Maybe
-- x = Just 42 -- значение внутри Maybe
-- применение:
-- maybeApply (Just fun) (Just x) = Just (fun x)
-- maybeApply _ _ = Nothing
--
-- print $ maybeApply f x  -> Just 43
-- print $ maybeApply Nothing x  -> Nothing
-- print $ maybeApply (Just (*2)) Nothing  -> Nothing
--
-- А для списков?
-- f = [(+1), (*2)]
-- x = [10, 20, 30]
-- что может означать listApply, когда есть список функций и список значений?
-- Применить каждую функцию к каждому элементу, получить список новых значений.
--
-- listApply fs xs = concatMap (\f -> map f xs) fs
-- print $ listApply f x  -> [11, 21, 31, 20, 40, 60]
-- Немного философии, почему для списков это логично делать именно так. Можно считать, что Maybe - это вычисления,
-- в которых иногда отсутствуют результаты. А в списках можно считать, что мы делаем вычисления, в которых может быть
-- несколько результатов.
-- В Haskell есть класс типов Аппликативные функторы:
-- class Functor f => Applicative f where
--   pure :: a -> f a -- вставить значение в функтор
--   (<*>) :: f (a -> b) -> f a -> f b -- та самая функция применения, рассмотренная выше
--   -- и ещё несколько, см. дальше
--
-- pure - это функция, которая помещает значение внутрь функции, "самым естественным способом". Что это значит будет
-- понятно из требований к аппликативному функтору, см. дальше.
-- Для Maybe: pure x = Just x, для списка pure x = [x]
--
-- Допустим, хотим сложить значения в двух Maybe. Или в двух списках сложить попарно значения. Как это будет написано?
-- x1 = Just 110
-- x2 = Just 654
-- надо сложить. Сначала давайте получим Just (+110), потом применим её к Just 654:
-- (+) <*> x1 <*> x2
--
-- Аппликативные функторы должны уметь делать не только <*>, но и функцию pure, которая вводит обычное значение в функтор.
-- 'pure 42' в типе Maybe это 'Just 42' . 'pure 42' в типе списка: [42]
--
-- Правила для функторов:
--   pure id <*> v = v
--   pure f <*> pure x = pure (f x)
--   и др.
--
--
-- #Монады
-- Вычисления с Maybe - это вычисления, в которых может не получиться результата. Если мы делаем последовательность
-- вычислений, в которой в какой-то момент результата не получилось, то и дальше считаем, что результата нет.
-- mSqrt :: Double -> Maybe Double
-- mSqrt x | x < 0 = Nothing
--         | otherwise = Just $ sqrt x
-- [хотим вычислить (sqrt(x) + 1) / 2]
-- eval x = let sx = mSqrt x
--              y = (+1) <$> sx
--            in (/2) <$> y
-- eval 9 -> Just 2.0
-- eval (-9) -> Nothing
--
-- sqrt(x) + sqrt(y) + 1
-- eval2 x y = let sx = mSqrt x
--                 sy = mSqrt y
--                 sq_plus_sq = Just (+) <*> sx <*> sy -- liftA2 (+) sx sy
--                in
--                 (+1) <$> sq_plus_sq
-- eval2 9 16 -> Just 8.0
-- eval2 (-9) 16 -> Nothing
--
-- sqrt(sqrt(x) - 3)
-- eval3 x = let sx = mSqrt x
--               sx_3 = (\t -> t - 3) <$> sx
--               -- ?? не получается сделать ни с <*>, ни с <$>
--               -- есть Just y и есть mSqrt :: Double -> Maybe Double
--               z = case sx_3 of
--                      Just n -> mSqrt n
--                      Nothing -> Nothing
--              in
--               z
-- eval3 100 -> Just ...
-- eval3 (-100) -> Nothing
-- eval3 4 -> Nothing
--
-- Хочется вспомогательную функцию типа <*> и <$>, чтобы она имела такой заголовок
-- bind :: Maybe a -> (a -> Maybe b) -> Maybe b
--
-- Если подумать, то
-- 1) fmap - частный случай этой функции. Сравним
--    bind :: Maybe a -> (a -> Maybe b) -> Maybe b
--    fmap :: Maybe a -> (a -> b) -> Maybe b
--    fmap f t = bind t (\x -> pure $ f x) или fmap f t = bind t (pure . f)
-- 2) <*> - тоже частный случай. Попробуйте её сам выразить через bind
--
-- Реализуем bind для Maybe:
-- maybeBind :: Maybe a -> (a -> Maybe b) -> Maybe b
-- maybeBind Nothing _ = Nothing
-- maybeBind (Just x) f = f x
--
-- снова sqrt(sqrt(x) - 3)
-- eval3 x = let sx = mSqrt(x)
--               sx_3 = maybeBind sx (\t -> Just (t - 3))
--              in
--               maybeBind sx_3 mSqrt
-- eval3 100 -> Just ...
-- eval3 (-100) -> Nothing
-- eval3 4 -> Nothing
--
-- короткий синтаксис
-- eval3 x = do
--             sx <- mSqrt x
--             sx_3 <- sx - 3
--             mSqrt sx_3
-- ответы те же
--
-- Давайте остановимся на одной вспомогательной функции
-- maybeBind :: Maybe a -> (a -> Maybe b) -> Maybe b
-- Функции даётся значение, очередной шаг вычисления, она возвращает новое значение.
-- maybeBind :: Maybe a -> (a -> Maybe b) -> Maybe b
-- maybeBind Nothing _ = Nothing
-- maybeBind (Just x) f = f x
--
-- взятие корня
-- mSqrt :: Double -> Maybe Double
-- mSqrt x | x < 0 = Nothing
--         | otherwise = Just $ sqrt x
-- sqrt(sqrt(x) - 3)
--
-- eval4 :: Double -> Maybe Double
-- maybeBind (Just x) mSqrt - это sqrt(x) внутри maybe
-- eval4 x = maybeBind (mSqrt x) (\a -> mSqrt (a - 3))
--
-- Эта программа похожа на последловательность шагов:
-- 1. mSqrt
-- 2. \a -> mSqrt (a - 3)
--
-- Посчитаем теперь sqrt x + sqrt y + 1
-- Шаги:
-- 1. sqrt x
-- 2. sqrt y
-- 3. sqrt x + sqrt y + 1
--
-- eval5 :: Double -> Double -> Maybe Double
-- eval5 x y = maybeBind (mSqrt x) (\a -> maybeBind (mSqrt y) (\b -> Just (a + b + 1)))
--
-- maybeBind - это функция '>>=', которая существует для Maybe и других типов.
-- eval5 :: Double -> Double -> Maybe Double
-- eval5 x y = (mSqrt x) >>= (\a -> (mSqrt y) >>= (\b -> Just (a + b + 1)))
--
-- для подообной записи есть do-нотация
-- eval5 :: Double -> Double -> Maybe Double
-- eval5 x y = do
--              a <- mSqrt x
--              b <- mSqrt y
--              return (a + b + 1)
-- Получился "императивный" код. return для Maybe это Just.
-- Перепишем теперь eval4
-- eval4 :: Double -> Maybe Double
-- eval4 x = (mSqrt x) >>= (\a -> mSqrt (a - 3))
-- или в do-нотации
-- eval4 x = do
--             a <- mSqrt x
--             mSqrt (a - 3) -- сразу возвращает Maybe
--
-- Комплексное число 1 :+ 2 это 1 + 2i
-- import Data.Complex
-- lSqrt :: Complex Double -> [Complex Double]
-- lSqrt x = [sqrt x, -sqrt x]
--
-- eval5 :: Complex Double -> Complex Double -> [Complex Double]
-- eval5 x y = do
--              a <- lSqrt x
--              b <- lSqrt y
--              return (a + b + 1)
--
-- Монада в хаскеле:
-- class Applicative a => Monad (m a) where
-- (>>=) :: m a -> (a -> m b) -> m b
-- return :: a -> m a -- превращает обычное значение в значение в монаде
-- (>>) :: m a -> m b -> m b
-- (>>) x y = x >>= (\_ -> y)
--
-- Пример монады. Вычисления с журналированием действий.
data Log a =
  Log a [String]
  deriving (Show) -- храним значение и список сообщений

-- Будем пользоваться функциями, которые при вычислениях могут выдать несколько сообщений:
add1 :: Int -> Log Int
add1 x = Log (x + 1) ["added 1"]

--
mul2 :: Int -> Log Int
mul2 x = Log (2 * x) ["multiplied by 2", "that was hard"]

--
-- Скажем, что Log это монада и объясним, как комбинировать вычисления.
instance Functor Log where
  fmap f (Log v msg) = Log (f v) msg

--
instance Applicative Log where
  pure x = Log x []
  Log f newMsg <*> Log v msg = Log (f v) (msg ++ newMsg)

--
instance Monad Log where
-- (>>=) :: Log a -> (a -> Log b) -> Log b
  (Log x messages) >>= f =
    let (Log y newMessages) = f x
      in Log y (messages ++ newMessages)

  -- return :: a -> Log a
  return x = Log x []

--
-- 2 + x + 1 через функции add1 и mul2
eval6 :: Int -> Log Int
eval6 x = do
  a <- mul2 x
  add1 a -- или b <- add1 a; return b
-- eval6 10 -- должен вернуть 21 и журнал ["умножила на 2", "добавила 1"]

--sumIO :: IO ()
--sumIO = do
--  a <- getLine
--  b <- getLine
--  print (read a + read b)
sumIO :: IO Int
sumIO = do
    a <- getLine
    b <- getLine
    return (read a + read b)


text:
{-# LANGUAGE TypeOperators #-}
module Types where

import           Data.Array.Repa
import           Data.Complex

type DIM6 = DIM5 :. Int

type R2S1Array   = Array U DIM3 (Complex Double) -- (Z :. numOrientations :. xLen :. yLen)
type R2S1RPArray = Array U DIM4 Double            -- (Z :. numOrientations :. numScales :. xLen :. yLen)
type R2S1T0Array = Array U DIM4 (Complex Double)  -- (Z :. numOrientations :. (L.length theta0freqs)  :. xLen :. yLen )
type R2Z1T0Array = Array U DIM4 (Complex Double) -- (Z :. (L.length thetafreqs) :. (L.length theta0Freqs) :. xLen :. yLen )

type R2Z1T0RadialArray = Array U DIM3 Double -- (Z :. (L.length thetafreqs) :. (L.length theta0Freqs) :. rLen)

type R2S1RPT0S0Array = Array U DIM6 (Complex Double)
-- (Z :. numOrientations :. numScales :. (L.length theta0freqs) :. (L.length scale0freqs)  :. xLen :. yLen )
type R2Z2T0S0Array = Array U DIM6 (Complex Double)
-- (Z :. (L.length thetafreqs)  :. (L.length scalefreqs) :. (L.length theta0Freqs) :. (L.length scale0Freqs) :. xLen :. yLen )

type R2Array   = Array U DIM2 (Complex Double) -- (Z :. xLen :. yLen)
type R2T0Array   = Array U DIM3 (Complex Double) -- (Z :. (L.length theta0freqs) :. xLen :. yLen)
type R2Z1Array   = Array U DIM3 (Complex Double) -- (Z :. (L.length thetafreqs) :. xLen :. yLen)
type R2Z2Array   = Array U DIM4 (Complex Double) -- (Z :. (L.length thetafreqs) :. (L.length scalefreqs) :. xLen :. yLen)

type R2T0S0Array   = Array U DIM4 (Complex Double) -- (Z :. (L.length theta0freqs) :. (L.length scale0freqs)  :. xLen :. yLen)

data R2S1RPPoint =
  R2S1RPPoint (Int, Int, Double, Double)
  deriving (Show, Read)

instance Eq R2S1RPPoint where
  (==) (R2S1RPPoint (x1, y1, _, _)) (R2S1RPPoint (x2, y2, _, _)) =
    x1 == x2 && y1 == y2

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field