1

in F#, how does one write a generic-math step function?

An (Oliver) Heaviside step function is function that returns zero if x is negative, otherwise it retuns one.

Here is a summary of my attempts so far:

// attempt 1:
let inline stepFct1< ^T when ^T : (static member op_GreaterThan: ^T * float -> bool) 
    >     (x:^T) : ^T = 
    //if (^T : (static member op_GreaterThan) (x 0.0) ) then x  //ouch fails also
    if  (>) x 0.0 then x
    else 0.0

compiler says: error FS0001: A type parameter is missing a constraint 'when ^T : comparison'

// attempt 2:
let inline stepFct2<^T when ^T : (static member (>): ^T * ^T -> bool) > (x:^T) : ^T = 
    match x with 
    | x when x > 0.0 -> 1.0
    | 0.0

FSC says: error FS0010: Unexpected infix operator in pattern

Motivation:

I am trying to rewrite Ian's Cumulative-Normal and Black-Scholes functions here to use Automatic Differentiation (DiffSharp). Ian's Cumulative Normal works on floats, I would like a generic version that works on any numeric type, including AutoDiff.DualG. The cumulative normal function contains a "greater than" statement.

EDIT: Gustavo, thanks, I have accepted your answer - the simple step function now compiles.

But it doesn't seem to help with the Cumulative Normal case. Given this code:

// Cumulative Normal Distribution Function - attempt to write a generic version
let inline CDF(x:^T) : ^T = 
    let (b1,b2,b3)  = (0.319381530, -0.356563782, 1.781477937)
    let (b4,b5)     = (-1.821255978, 1.330274429)
    let (p , c )    = (0.2316419  ,  0.39894228)
    let (zero, one) = (LanguagePrimitives.GenericZero, LanguagePrimitives.GenericOne)
    if x > zero then
        let t = one / (one + p * x) 
        (one - c * exp( -x * x / 2.0)* t * (t*(t*(t*(t*b5+b4)+b3)+b2)+b1)) 
    else
        let t = 1.0 / (one - p * x) 
        (c * exp( -x * x / 2.0)* t * (t*(t*(t*(t*b5+b4)+b3)+b2)+b1))

FSI says:

C:\stdin(116,32): warning FS0064: This construct causes code to be less generic 
than indicated by the type annotations. 
The type variable 'T has been constrained to be type 'float'.

val inline CDF : x:float -> float
> CDF 0.1M;;
CDF 0.1M;;
----^^^^
C:\stdin(122,5): error FS0001: This expression was expected to have type
    float
but here has type
    decimal
>

Does anyone know how to make CDF generic?

Richard
  • 106,783
  • 21
  • 203
  • 265
Piga-fetta
  • 119
  • 1
  • 6
  • 1
    Your update is really a new question, but you have missed the point of the answer. If you have a floating point literal anywhere, the code is no longer generic. You will need to find generic ways to construct the constants. – John Palmer Nov 28 '14 at 04:26
  • @JohnPalmer is right, that's more challenging. The lib I mentioned provides a way by using rational numbers as input and converting them to the destination type. See the update in the answer. – Gus Nov 28 '14 at 09:51

2 Answers2

6

Use LanguagePrimitives.GenericZero / GenericOne and let type inference do the rest

// attempt 1:
let inline stepFct1 x =
    let zero = LanguagePrimitives.GenericZero 
    if x > zero then x
    else zero

I had a look at the link you sent with the function you want to implement. FSharpPlus (F#+) may help you to write generic math code since it contains a dedicated Generic Numbers module. Or at least you can grab some techniques from there.

UPDATE

Regarding your updated question, which takes the complexity to a higher level, here is a solution using the latest version of F#+:

let inline CDF(x: ^T) : ^T = 
    let num x = fromRational (x </ratio/> 1000000000I)
    let (b1,b2,b3)  = (num 319381530I   , num -356563782I  , num 1781477937I)
    let (b4,b5)     = (num -1821255978I , num 1330274429I)
    let (p , c )    = (num  0231641900I , num 0398942280I)
    let (zero, one, two) = 0G, 1G, 2G
    if x > zero then
        let t = one / (one + p * x) 
        (one - c * exp( -x * x / two)* t * (t*(t*(t*(t*b5+b4)+b3)+b2)+b1)) 
    else
        let t = one / (one - p * x) 
        (c * exp( -x * x / two)* t * (t*(t*(t*(t*b5+b4)+b3)+b2)+b1))

But fromRational and Ratio need to be defined, I did a working example here so you can test your function, which works nicely with float and float32.

If you are interested in Generic Maths feel free to contribute with code or use cases.

Gus
  • 25,839
  • 2
  • 51
  • 76
  • see the 'EDIT' above, Gustavo's approach works for a simple step function. Answer accepted, but does not allow writing of a generic CDF. – Piga-fetta Nov 28 '14 at 03:47
  • Also thanks for the GenericMath link above. This is the first generic-math library I see - I will try to compile it tomorrow – Piga-fetta Nov 28 '14 at 04:07
  • Hi Gustavo - I just ran your version of CumulativeNormal and I have two comments: 1) it is hideous - all float constants replaced by huge integers. 2) it is amazing. Until now, all the numeric libraries I have seen are either targeting a single type (double), or they are implemented as C++ headers. Your idea allows for a generic library - could be revolutionary! – Piga-fetta Nov 28 '14 at 21:28
  • @Piga-fetta Thanks. The inputs are expressed as two bigIntegers which represents a rational, that way there is no precision loss in the input, but then they are immediately converted to the target type, before the calculations. I'm happy to know that this GenericMath module can be useful for you, I will spend more time on it, but I will need some sample code. If you have a Github account, let's continue the discussion in the project page, your feedback will be much appreciated. – Gus Nov 28 '14 at 21:57
0

With addition of Generic Math to .NET 7 and latest F# you can do something like this (I was not able to remove and 'T:comparison constraint though comparison operators are defined for INumber, submitted bug for the compiler):

let inline stepFct<'T when 'T :> INumber<'T> and 'T:comparison> x:'T = 
    if x > 'T.Zero then x
    else 'T.Zero

Similar approach can be used for second one:

let inline stepFct<'T when 'T :> IExponentialFunctions<'T> and  'T : comparison> (x:^T) : 'T=
    let (b1,b2,b3)  = ('T.CreateChecked(0.319381530), 'T.CreateChecked(-0.356563782), 'T.CreateChecked(1.781477937))
    let (b4,b5)     = ('T.CreateChecked(-1.821255978), 'T.CreateChecked(1.330274429))
    let (p , c )    = ('T.CreateChecked(0.2316419)  ,  'T.CreateChecked(0.39894228))
    let (zero, one) = ('T.Zero, 'T.One)
    if x > zero then
        let t = one / (one + p * x) 
        (one - c * exp( -x * x / 'T.CreateChecked(2.0))* t * (t*(t*(t*(t*b5+b4)+b3)+b2)+b1)) 
    else
        let t = one / (one - p * x) 
        (c * exp( -x * x / 'T.CreateChecked(2.0))* t * (t*(t*(t*(t*b5+b4)+b3)+b2)+b1))
Guru Stron
  • 102,774
  • 10
  • 95
  • 132