Please give feedback on my function

I’m very new to Julia, and to get started I’m translating some old MATLAB functions into a package that will be used analyze glacier measurements. Here’s my first attempt at writing a function, and I want to make sure I’m starting out on the right foot.

  • Does anything need to be changed to make it more Julia-y?
  • Does anything–even the docstring formatting–jump out at you as bad syntax?

Feel free to be brutal in your feedback, but please let me know which suggestions are essential versus which might be minor personal preferences.

"""
    T = freeboard2thickness(F,kwargs)
  
Estimate ice thickness from surface height above sea level, assuming hydrostatic equilibrium.
  
# Input
- `F` (m) freeboard surface height above sea level.
- `rho_ice` (kg/m³) density of ice layer  (default: `917`).
- `rho_seawater` (kg/m³) density of seawater   (default: `1027`).
- `rho_snow` (kg/m³) density of snow layer (default: `350`). Inconsequential unless `T_snow>0`.
- `T_snow` (m) thickness of snow layer   (default: `0`).
  
# Output
- `T` (m) thickness of ice layer. 
  
# Examples
### Example 1: Simple case
The surface of an iceberg is 80 cm above the surrounding ocean surface, and the iceberg is 
made of pure ice, without any snow on it. How thick is the iceberg? 
```julia 
julia> freeboard2thickness(0.8)
7.46909090909091
```
### Example 2: Specifying parameters 
The surface of an iceberg is 80 cm above sea level, and has a 40 cm layer of snow on top, 
whose density is 300 kg/m³. How thick is the ice portion of this iceberg? 
```julia
julia> freeboard2thickness(0.8; T_snow=0.4, rho_snow=300)
4.825454545454546
```
# Author Info 
Chad A. Greene, NASA Jet Propulsion Laboratory 
April, 2022
"""
```julia
function freeboard2thickness(
    F; 
    rho_ice=917, 
    rho_seawater=1027, 
    rho_snow=350,
    T_snow=0
    )
    T = F*rho_seawater/(rho_seawater-rho_ice) - T_snow*(rho_seawater-rho_snow)/(rho_seawater-rho_ice); 
    return T
end
```
1 Like

Covered by this PSA.

(Edit: also changed the title, don’t know if insult is the easier way to get people interested. :sweat_smile:)

2 Likes

As for your function, it is all right. :slight_smile:

If you want and there is an established well known nomenclature on your area, you can use Unicode for the \rho's and alike, some people like it, some don’t. One thing to favor not using Unicode is that not everybody has a code viewer/font that renders well Unicode, so

For the docs

don’t just put kwargs, you can use the way it is used in Base, this is

    freeboard2thickness(freeboard; rho_ice=917, rho_seawater=1027, rho_snow=350, T_snow=0)

and later on explain from where you got those default parameters. :slight_smile:

1 Like

This is what I would write for docstrings + function, I did the next changes

  • Changed all \rho for density, this way you don’t need to specify the meaning, as it is clear by the name
  • change the name of the function, it is more readable for non-native English speakers (I’m one :sweat_smile:)
  • changed keyword arguments to be similar to base syntax of docstrings (you can see ?sort)
  • reduce the examples, as I don’t think they need to be so verbose (not a critic, just my opinion), specially because the keyword arguments are more expressive

Personally, I wouldn’t mind putting the author there, as that information would come in the package registry and license. But does not harm in any way. :slight_smile:

For the code, I just separate it in two lines for it not to be too large. The code itself is alright, I just will mention that maybe you want to use Unitful.jl for the units. :smiley:

"""
        freeboar_to_thickness(freeboard_thickness; ice_density=917, seawater_density=1027, snow_density=350, thickness_snow=0)

Estimate ice thickness (in meters) from surface height above sea level, assuming hydrostatic equilibrium. Thickness are expected to be in meters and densities in kg/m^3.

The default parameters come from (here you put the reference).

# Examples

```
julia> freeboard_to_thickness(0.8) # pure ice, 80cm above sealevel
7.46909090909091

julia> freeboard_to_thickness(0.8; snow_thickness=0.4, snow_density=300) # 80cm above sea level, 40cm of snow with density 300 kg/m^3
4.825454545454546
```

# Author Info

Chad A. Greene, NASA Jet Propulsion Laboratory
April, 2022
"""
function freeboard_to_thickness(
    freeboard; 
    ice_density=917, 
    seawater_density6=1027, 
    ρsnow=350,
    thickness_snow=0
    )
    thickness = freeboard * seawater_density / (seawater_density - ice_density) - # 
        snow_thickness * (seawater_thickness - rho_snow) / (rho_seawater - rho_ice); 
    return thickness
end

This is what it looks like in my terminal

1 Like

A couple more suggestions:

  • The semicolon ending the line T = ...; is unnecessary and unjulian.
  • The line containing return is unnecessary, since in its absence the function will return the value of the final expression (T in this case). It’s not wrong, but after writing lots of Julia functions, most people don’t bother with a return statement if it is the final statement of the function.
2 Likes

It may be common practice, but the Blue Style style guide says to always explicitly return a value for longform method definitions. If the method has no return value, it should still explicitly return nothing.

Implicit returns can be confusing to read, so I would agree that it’s best to always return something

9 Likes

Note, OP, that this requires no change in the definition of the function. If you call this function with compatible unitful quantities, the result will be a unitful quantity too.

3 Likes

But the default parameters would need to be uniful or you’d have to always explicitly give them, no?

Correct. Since you cannot subtract a dimensionless number from a density, for instance.

1 Like

Good to know. So when should I be using a semicolon?

I wouldn’t define thickness there, just return the value directly. We already know it’s a thickness from the name of the function and docstring.

1 Like

There might be a way to calculate the correct units for the default parameter values with Unitful.jl in terms of the positional argument’s type, too, but I don’t know the subject matter to do that.

I’m not sure I follow. Can you explain this point a little more?

Also, I see you removed the Input and Output sections of the docstring in your suggested rewrite. Why is that?

Semicolons are typically only used to terminate lines in the REPL or in a notebook when you wish to suppress printout of the expression’s value. These printouts don’t occur for code executed from a file.

3 Likes

Sure,

f(x; a=3) = a*x

f(x::Unitful.Quantity; a=3u"m") = a*x

The given suggestions are great so far.
I would like to make another one that would make your code slightly more julian.

If I inderstand your function correctly, what you compute is the thickness of the glacier from the freeboard?
Are there other things you can get glacier thicknesses from? If the answer is yes, I would do the following:

function glacier_thickness(f:: Freeboard)
    F =  f.F
    rho_ice = f.rho_ice
    rho_seawater =f.rho_seawater
    rho_snow = f.rho_snow
    T_snow = f.T_snow

    return F*rho_seawater/(rho_seawater-rho_ice) -    T_snow*(rho_seawater-rho_snow)/(rho_seawater-rho_ice); 
end

Where Freeboard is a struct that contains all data:

Base.@kwdef struct Freeboard
    F
    rho_ice=917 
    rho_seawater=1027
    rho_snow=350
    T_snow=0
end
   

Maybe in the function you can use some macro to unpack the contents of the struct to avoid all the boilerplate.
Anyway, what this does is that now you can potentially make a second method (in julia parlance) of glacier_thickness that acts on a different object like for a glacier on land for instance?
Like this you leverage what’s called multiple dispatch.
This is the way to make generic code since your function glacier_thickness can be used on any object of relevance. You don’t need anylonger to define a function continentalfreeboard2thickness.
You could just do

function glacier_thickness(c::ContinentalFreeboard)
#  I do stuff here
end

I mean something more like

julia> f(x; y=3*unit(x)) = x + y
f (generic function with 1 method)

julia> f(10u"cm")
13 cm

I removed output and input because I think is clear by the arguments what the input is, and by the description what the output is. This makes the docstrings more concrete and shorter.

By saying that I changed the styles of kwargs to match Base, I was meaning that the default values are specified directly in the function signature, this is how it is done in the docs for sort, as you may see in the REPL.

This also makes for the input section to be redundant, as the signature says the default values and the meaning of kwargs obvious. I know this is concise but not the clearer (for a beginner), but I do think this convention is better for the REPL help/docstring and something more elaborated to be in the documentation. :slight_smile:

1 Like

As of version 1.7, Julia has syntax for unpacking object properties:

julia> struct A
           x
           y
       end

julia> a = A(1, 2)
A(1, 2)

julia> (; x, y) = a
A(1, 2)

julia> x
1

julia> y
2

(Also works for named tuples, or any object with a getproperty method.)

That only works if x and y necessarily has the same unit, and kind of removes the point of dimensional analysis. Say we have f(x; y=3) = x + y, and we want to change y to default to 3 m, then we want f(1u"s") to fail, and f(1u"cm") to give 3.01 m.