13
$\begingroup$

I have a function value like this:

groups = N[{{40, 23, 213, 79, 18, 152, 105, 191, 119, 243}, {86, 3, 152, 89, 
     52, 36, 132, 196, 159, 108, 200}, {253, 223, 167, 197, 147, 112, 
     225, 98, 159, 163, 245, 240, 267}, {89, 164, 156, 124, 15, 253, 
     41, 171, 209, 132}}];

PopulationVariance1[list_] := Mean[(list - Mean[list])^2]

value[groups_] := Module[{mean, SSA, SSE}, mean = Mean[Flatten[groups]];
  SSA = Total[Table[Length[group]*(Mean[group] - mean)^2, {group, groups}]];
  SSE = Total[Table[Length[group]*PopulationVariance1[group], {group, groups}]];
  (SSA/(Length[groups] - 1))/(SSE/(Length[Flatten[groups]] - Length[groups]))]

value[groups]

3.6118

This is my current compilation try:

dec = FunctionDeclaration[PopulationVariance, 
   Typed[{"PackedArray"::["Real64", 1]} -> "Real64"]@
    Function[list, Mean[(list - Mean[list])^2]]];

cfvalue = FunctionCompile[dec, 
  Function[Typed[groups, "ListVector"::["ListVector"::["Real64"]]], 
   Module[{mean, SSA, SSE}, 
    mean = Mean[
      Cast[Catenate[groups], "PackedArray"::["Real64", 1]]];
    SSA = Total[Cast[
       Length[#]*(Mean[Cast[#, "PackedArray"::["Real64", 1]]] - 
             mean)^2 & /@ groups, "PackedArray"::["Real64", 1]]];
    SSE = Total[Cast[
       Length[#]*PopulationVariance[
           Cast[#, "PackedArray"::["Real64", 1]]] & /@ groups, 
       "PackedArray"::["Real64", 1]]];
    (SSA/(Length[groups] - 1))/(SSE/(Length[Catenate[groups]] - 
         Length[groups]))]]]

cfvalue[groups]

3.6118

It works indeed. However, I found that the compiler function was even slower than the uncompiled one.

lis = DeleteCases[Table[ResourceFunction["RandomSplit"][RandomReal[255, 1000], 10], 
    100000], {}, {2}];
Timing[value /@ lis;](*{15.2813, Null}*)
Timing[cfvalue /@ lis;](*{23., Null}*)

Could you please help me improve the speed of the cfvalue?


$Version

13.1.0 for Microsoft Windows (64-bit) (June 16, 2022)

$\endgroup$
7
  • $\begingroup$ Dear @yode. I get syntax error for your Typed[{"PackedArray"::["Real64",1]}->"Real64"] in Version 12.3. Concerning your question, I see that your value code is nice and simple. I wonder if in 202x, with just-in-time compilation and all that jazz, Mathematica does not automatically optimize code such as yours, and nothing is gained by explicitly compiling yourself? I am sure you know much more about this than I do. $\endgroup$
    – user293787
    Jul 20, 2022 at 12:51
  • 3
    $\begingroup$ @user293787 The new syntax of TypeSpecifier is added in v13.1. $\endgroup$
    – xzczd
    Jul 20, 2022 at 12:56
  • 1
    $\begingroup$ @user293787 Hi, I'm in version 13.1,If you are in 12.3, you should use Typed[{TypeSpecifier["PackedArray"]["Real64", 1]} -> "Real64"] $\endgroup$
    – yode
    Jul 20, 2022 at 13:10
  • $\begingroup$ On my computer, the whole operation with value on the lis took 19 seconds (cfvalue was 32 seconds), I was curious to see how a low-level version would perform, so I tried it on Rust and it took 5.6 seconds which if we ignore the time it takes to convert the input to NumericArray, it was 0.4 second! Clearly, there is much room for improvement. $\endgroup$
    – Ben Izd
    Jul 21, 2022 at 5:22
  • $\begingroup$ @BenIzd Thanks, this confirms that OP has a legit expectation that the compiled code should be faster. I would find it useful if you posted details about this as an answer, even if it is more supporting material rather than an actual answer. But not sure if Mathematica SE allows such answers. $\endgroup$
    – user293787
    Jul 21, 2022 at 5:40

3 Answers 3

13
$\begingroup$

The Code

Here's some equivalent code that performs a better:

cfValue = 
 FunctionCompile[
  Function[Typed[groups, "ListVector"::["PackedArray"::["Real64", 1]]],
   Module[{
     PopulationVariance = 
      Function[{Typed[g, "PackedArray"::["Real64", 1]]},
       Mean[(g - Mean[g])^2]
       ], mean, total, SSA, SSE, len},
    
    total = 
     Fold[Plus[#1, #2] &, 
      Map[Typed[Total, {"PackedArray"::["Real64", 1]} -> "Real64"], 
       groups]];
    len = 
     Fold[Plus[#1, #2] &, 
      Map[Typed[
        Length, {"PackedArray"::["Real64", 1]} -> "Integer64"], 
       groups]];
    
    mean = total/ len;
    
    SSA = 
     Fold[Plus[#1, #2] &, Length[#]*(Mean[#] - mean)^2 & /@ groups];
    SSE = 
     Fold[Plus[#1, #2] &, 
      Map[Length[#]*PopulationVariance[#] &, groups]];
    (SSA/(Length[groups] - 1))/(SSE/(len - Length[groups])) 
    ]
   ]
  ]

Timing

On my machine with Mathematica 13.1:

In[111]:= AbsoluteTiming[cfValue[#] & /@ lis;]

Out[111]= {3.04551, Null}

In[112]:= AbsoluteTiming[value[#] & /@ lis;]

Out[112]= {10.4156, Null}

Why

The casting in code from the question was cloning the data each time which is very expensive.

This function avoids casting by instead using a Fold with Plus to total the ListVectors.

$\endgroup$
12
  • 2
    $\begingroup$ Thank you for the clarification! (+1) A bit off-topic: does casting in the new compiler always cause data copy (even the "BitCast" version)? $\endgroup$
    – Silvia
    Jul 21, 2022 at 10:48
  • 2
    $\begingroup$ It depends upon the types and the scenario. A cast of an atomic type to another atomic type (eg, Integer32 to Integer64) necessarily will. A cast from PackedArray[Int32] to NumericArray[Int32] it is possible the data would not be cloned as they are both immutable arrays. However, as PackedArray and ListVector are fundamentally different data types (one an array, one a linked list or something) there is necessarily a clone of the data. This is the same as if you cast a list to an array in any language. $\endgroup$
    – A.Richards
    Jul 21, 2022 at 11:26
  • 2
    $\begingroup$ Hi Yode, Yes I added the #1 and & when I was struggling with the type system but its unnecessary, whoops. The Total does not work because Total has not been implemented for ListVectors. Hopefully it will be added for a future update. $\endgroup$
    – A.Richards
    Jul 21, 2022 at 15:26
  • 2
    $\begingroup$ @A.Richards Thank you for that very useful information! @yode You can use $CompilerEnvironment["TypeEnvironment"]["getPolymorphicList", f ] to check compilable type of function f. It shows Total is only defined for container types belong to DenseArrays, OTOH Fold is defined for all types belong to RangeIterable. We can check ListVector is not a kind of DenseArrays but is a kind of RangeIterable. So Fold is compilable here but Total is not. $\endgroup$
    – Silvia
    Jul 22, 2022 at 8:50
  • 2
    $\begingroup$ @yode BTW to check if one type is a "kind" of certain abstract type, you can do this: env=$CompilerEnvironment@"TypeEnvironment";env["implementsQ","resolve"~env~"ListVector","resolve"~env~TypeFramework`AbstractType@"RangeIterable"] . $\endgroup$
    – Silvia
    Jul 22, 2022 at 8:55
9
$\begingroup$

I decide to use Rust to show how much improvement we can achieve while also benefiting from the Rust memory management system which prevents us to leak memory (to some extent).

Recently Wolfram started a project called wolfram-library-link-rs which is a Rust crate (aka library/module/etc) which allows users to code in Rust and export their functions without doing much for the interface with Mathematica.

To follow this workflow you need to install Rust and have an internet connection when you build the project (I think only the first time is required).

Create a Rust project, replace Cargo.toml content with the following:

[package]
name = "population_variance"
version = "0.1.0"
edition = "2021"

[lib]
crate-type = ["cdylib"]

[dependencies]
wolfram-library-link = "0.2.5"

and replace src/lib.rs content with this:

use wolfram_library_link::{export, NumericArray};

type Precision = f64;
type ArrayType = Precision;

fn population_variance(array: &[ArrayType]) -> Precision {
    let mean = array.iter().sum::<ArrayType>() as f64 / array.len() as Precision;
    array.iter().map(|&x| (x as Precision - mean).powf(2.0)).sum::<Precision>() as Precision / array.len() as Precision
}

fn value(groups: Vec<&[ArrayType]>) -> Precision {
    let total: ArrayType = groups.iter().map(|group| group.iter().sum::<ArrayType>()).sum::<ArrayType>();
    let flatten_length: usize = groups.iter().map(|group| group.len()).sum();
    let mean = total as Precision / flatten_length as Precision;
    let groups_len = groups.len();

    let ssa = groups.iter().map(|group| group.len() as Precision * ((group.iter().sum::<ArrayType>() as Precision) / (group.len() as Precision) - mean).powf(2.0)).sum::<Precision>();
    let sse = groups.iter().map(|group| group.len() as Precision * population_variance(group)).sum::<Precision>();
    (ssa / (groups_len - 1) as Precision) / (sse / (flatten_length - groups_len) as Precision)
}

#[export]
fn mathematica_interface(data: &NumericArray<ArrayType>, lengths: &NumericArray<i64>) -> Precision {
    let data_slice = data.as_slice();
    let mut vectors: Vec<&[ArrayType]> = Vec::with_capacity(lengths.flattened_length());
    let mut current_index: usize = 0;
    for &length in lengths.as_slice() {
        let length = length as usize;

        vectors.push(&data_slice[current_index..(current_index + length)]);
        current_index += length;
    }

    value(vectors)
}

When you build the project (with build --release), you will get a population_variance.dll in the target/release directory.

In Mathematica use the following code to use the library (change the path to your address):

compiledValue = 
 LibraryFunctionLoad["C:\\population_variance.dll", 
  "mathematica_interface", {{LibraryDataType[NumericArray, "Real64"], 
    "Constant"}, {LibraryDataType[NumericArray, "Integer64"], 
    "Constant"}}, Real];

ClearAll[RustValue];

RustValue[data_] := Block[{lengths = Length /@ data},
  compiledValue[NumericArray[Flatten@data, "Real64"], 
   NumericArray[lengths, "Integer64"]]
  ]

Now it will work as a normal function:

Block[{data = N@{{1, 2, 3}, {4, 5, 6, 7, 9}}},
 {value[data], cfvalue[data], RustValue[data]}
 ]

(* Out: {11.8125, 11.8125, 11.8125} *)

Testing larger dataset:

Block[{data = N@lis},
 {value /@ data // Total // AbsoluteTiming, 
  cfvalue /@ data // Total // AbsoluteTiming, 
  RustValue /@ data // Total // AbsoluteTiming}
 ]

(* Out: {{19.3707, 100129.}, {32.3437, 100129.}, {5.59327, 100129.}} *)

If you want to measure only the calculation on compiled code without including the time needed for converting the input:

Block[{a = NumericArray[Flatten@#, "Real64"] & /@ lis, 
   b = NumericArray[Length /@ #, "Integer64"] & /@ lis},
  MapThread[compiledValue, {a, b}] // Total // AbsoluteTiming
  ] // AbsoluteTiming

(* Out: {5.96117, {0.416077, 100129.}} *)

The first timing (5.96) is for the whole operation, the second (0.41) is for the compiled code, and 100129. exist to show it returns the same result as shown before.

How did we handle ragged arrays

As you might have seen, our compiled function has two arguments, first the flattened data and second the length of each list. In Rust, we partitioned the data like TakeList.

Note that I still have so many things to learn and I'm sure the code above can be further improved.

$\endgroup$
2
  • $\begingroup$ Is my c++ more efficient than your rust? $\endgroup$
    – yode
    Jul 31, 2022 at 16:03
  • $\begingroup$ @yode Yes, mine took 5.6 seconds which I could reduce to 4.8 but yours is 4.3 seconds on my computer. I guess I should start learning Chinese ;) $\endgroup$
    – Ben Izd
    Aug 1, 2022 at 4:20
3
$\begingroup$

Add a method from this post:

Needs["CCompilerDriver`"]
lib=CCompilerDriver`CreateLibrary["
#include<iostream>
#include<vector>
#include<algorithm>
#include<numeric>
#include \"WolframCompileLibrary.h\"

using namespace std;

struct ManagedArray{
    MTensor* data;
    mint refcount;
};

struct FixedArray{
    ManagedArray data;
    mint refcount;
};

struct RaggedList{
    FixedArray* data;
    mint refcount;
};

mreal value(vector<vector<mreal>> pix_values) {
    vector<mreal> flatten_pix_values;
    for (auto const& group : pix_values) {
        flatten_pix_values.insert(flatten_pix_values.end(), group.begin(), group.end());
    }

    mreal mean_value = std::accumulate(begin(flatten_pix_values), end(flatten_pix_values), 0.0) / flatten_pix_values.size();

    mreal SST = 0.0;
    for_each(begin(flatten_pix_values), end(flatten_pix_values), [&SST, mean_value](const mreal d) {SST += pow(d - mean_value, 2); });

    mreal SSA = 0;
    for (auto const& group : pix_values) {
        SSA += group.size() * pow(std::accumulate(begin(group), end(group), 0.0) / group.size() - mean_value, 2);
    }
    mreal SSE = SST - SSA;
    int a = pix_values.size() - 1, b = flatten_pix_values.size() - pix_values.size();

    return mreal((SSA / a) / (SSE / b));
}

EXTERN_C DLLEXPORT mreal ragged_value(RaggedList inlis) {
    vector<vector<mreal>> vecs;
    
    for (mint i = 0; i < inlis.data->data.refcount; i++) {
        MTensor tensor = inlis.data->data.data[i];
        auto ptensor = (mreal*)tensor->data;
        vector<mreal> vec(ptensor, ptensor + tensor->nelems);
        vecs.push_back(vec);
    }

    return value(vecs);
}","Ragged_Value",Language->"C++"];
dec=LibraryFunctionDeclaration[fun->"ragged_value",lib,{"ListVector"::["PackedArray"::["Real64",1]]}->"Real64"];
fcvalue=FunctionCompile[dec,Function[{Typed[lis,"ListVector"::["PackedArray"::["Real64",1]]]},fun[lis]]]

Performance

lis = DeleteCases[
   Table[ResourceFunction["RandomSplit"][RandomReal[255, 1000], 10], 
    100000], {}, {2}];
RepeatedTiming[value /@ lis;]
RepeatedTiming[fcvalue /@ lis;]

{11.0433, Null}

{2.42399, Null}

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.