text:
Environ Eng Res > Volume 25(3); 2020 > Article
Adeogun, Ganiyu, Ladokun, and Ibitoye: Evaluation of hydrokinetic energy potentials of selected rivers in Kwara State, Nigeria

### Abstract

This Hydrokinetic energy system is the process of extracting energy from rivers, canals and others sources to generate small scale electrical energy for decentralized usage. This study investigates the application of Soil and Water Assessment Tool (SWAT) in Geographical Information System (GIS) environment to evaluate the theoretical hydrokinetic energy potentials of selected Rivers (Asa, Awun and Oyun) all in Asa watershed, Kwara state, Nigeria. SWAT was interfaced with an open source GIS system to predict the flow and other hydrological parameters of the sub-basins. The model was calibrated and validated using observed stream flow data. Calibrated flow results were used in conjunction with other parameters to compute the theoretical hydrokinetic energy potentials of the Rivers. Results showed a good correlation between the observed flow and the simulated flow, indicated by ash Sutcliffe Efficiency (NSE) and R2 of 0.76 and 0.85, respectively for calibration period, and NSE and R2 of 0.70 and 0.74, respectively for the validation period. Also, it was observed that highest potential of 154.82 MW was obtained along River Awun while the lowest potential of 41.63 MW was obtained along River Asa. The energy potentials obtained could be harnessed and deployed to the communities around the watershed for their energy needs.

### 1. Introduction

Hydrokinetic is the process of extracting energy from ocean currents, waves, rivers, canals, tidal ocean and other water bodies to generate electric power. Hydrokinetic and wave energy technologies require no impoundments, so that the challenges and externalities associated with dam construction and reservoir creation are not relevant in its implementation [1]. The hydrokinetic energy conversion system has been described as the type of system that converts the kinetic energy of flowing water bodies with little or no head to generate other useable form of energy [2]. The available hydrokinetic power potentials that can be obtained from a river depend on the speed of the river, ocean, or tidal current and this is a function of the density of water and speed of the current cubed [3].
The breakthrough in this technology suggested that kinetic energy of flowing water can be used to generate electricity in flowing water without stress. Hydrokinetic energy is one of the newest clean energy that is derivable from the power of moving water and is unhindered by many of the weaknesses experienced with other clean energy sources [4]. Hence, globally, efforts have been channeled in this direction towards harnessing potential of rivers using hydrokinetic energy technology [57]. Natural Resources Canada [8] has completed a study to develop a methodology to assess the hydraulic kinetic energy contained in Canadian rivers. Also, in Nigeria, few studies have been initiated in this direction towards its development and implementation, for example, [9] conducted a field study at a selected site along Ero River, North-central Nigeria to assess the hydrokinetic energy conversion potentials of the river reach. Variants of data including hydrology, topography, bathymetry and digital elevation models (DEM) were used in conjunction with year round on-site measurements for the determination of the hydraulic and hydrological parameters. These were used to obtain the instantaneous in-stream velocities necessary for the computation of the hydrokinetic power potential at selected location.
The application of soil and water assessment tool (SWAT) in hydrological analysis of upper catchment of Ikere Gorge Basin for hydrokinetic energy estimation was examined by [10]. The total theoretical hydrokinetic energy potential of the 10 basin selected was estimated as 36.4 MW. Also, potential hydrokinetic was estimated to have highest value in Oshe at Onikankan (9.542 MW) and lowest in Kojuoba at Olonje (0 MW).
The assessment of the Hydropower Potential of Kangimi Reservoir (KR) in Kaduna State Nigeria was conducted by [11]. The study evaluated the hydropower potentials of the KR and the result showed that KR can potentially generate 1.109 MW and 0.692 MW if the KR dam axis is, respectively placed at 612 m and 604 m above mean sea level (AMSL).
Despite the fact that Nigeria is one of the most populous countries in Africa with a population of over 180 million [12], the informal sector of the economy that accommodates small scale industries that cater for majority of the people heavily depends on epileptic supply of power for economic production. In spite of the huge amount of capital invested in the generation and distribution of electrical energy in Nigeria on a yearly basis, quite a lot of people especially, the rural dwellers, most of which are the custodian of the cottage industries, are not privileged to be connected to the grid or have direct access to electricity. Additional problems such as vandalization of gas pipe line, seasonal wind variation, and obsolete machinery in the hydropower dams resulted in uneven distribution of power and sabotage of various sources of power generation in Nigeria [13].
Therefore, alternative energy generation to augment the generating capacities of existing hydropower plants and thermal stations to have stable power supply in Nigeria is inevitable. It is also a matter of urgency for energy stakeholders in Nigeria to harness the existing potentials of water resources in line with the changing world of hydrokinetic technology in the power generation as obtainable in many part of the world. A shift from other sources of energy to recent technology of hydrokinetic power production will add more values to economic growth and development of both the private and public sectors. According to Federal Ministry of Water Resources [14], small hydropower schemes have been integrated into some dam projects across the country to increase energy supply of the nation. However, the hydrokinetic potentials of river currents in Nigeria have been, to date, largely unexplored [15].
Consequently, this research examined the potentialities of using an hydrological model, SWAT embedded in GIS environment to model the hydrology of Asa watershed in Kwara State, Nigeria and used the modeling results in conjunction with other parameters to estimate the theoretical hydrokinetic energy potential of selected Rivers in the catchment.

### 2.1. Study Area

The study area is located between Latitude 8036′ and 8024′ North and Longitudes 4036′ and 4010′ East with estimated total area of about 639,239.07 Ha and it lies within Kwara State in Nigeria. The topography of the area could be described as a fair representative of surrounding plains with undulating and very broad and gentle slopes. The study area is characterized with an altitude between 265 m and 457 m AMSL.
The climate of the watershed could be described as humid tropic type and characterized by both wet and dry seasons with a mean annual temperature that ranges from 25–28.9°C during the wet season and 33–37°C in the dry season [16]. The watershed has an annual mean rainfall of about 1,150 mm. Some of the major rivers within the catchment include River Asa, River Oyun, and River Awun. Fig. 1 shows the location Map of the study area.

### 2.2. SWAT Model Selection and Description

In this study, SWAT was selected as the modeling tool. The tool was developed by the United States Department of Agriculture to predict the impact of land management practices on water, sediment and agricultural chemical yields in large un-gauged basins [17]. The modeling tool was interfaced with Mapwindow GIS to simulate the hydrology and predict the flow into sub-basins in the watershed. The SWAT model is a long-term, continuous simulation watershed model which operates on a daily time step and is designed to predict the impact of management on water and sediment. The model is physically based, computationally efficient, and capable of simulating a high level of spatial details by allowing the division of watersheds into smaller sub watersheds.
SWAT incorporates the effects of weather, surface runoff, evapo-transpiration, irrigation, sediment transport, groundwater flow, crop growth, nutrient yielding, pesticide yielding and water routing, as well as the long term effects of varying agricultural management practices [18, 19]. In the SWAT model, the watershed is partitioned into sub basins that are further subdivided into one or several homogeneous hydrological response units (HRUs) with relatively unique combinations of land cover, soil and topographic conditions. The hydrological component of the model calculates soil-water balance at each time step based on daily amounts of precipitation, runoff, evapo-transpiration, percolation and base flow.

### 2.3. Model Input Data Collection and Processing

The spatially distributed data (GIS input) needed for the running of Mapwindow SWAT include the DEM, soil data, land use and stream network layers. Weather data are also necessary for prediction of stream flow in the watershed. The GIS interface of the model was used to discretize the catchment area and geo-processed the spatial data into formats required by the model.

#### 2.3.1. DEM

Topography of a place is defined by a DEM that describes the elevation of any point in a given area at a specific spatial resolution. 90 m by 90 m resolution DEM of the study area was obtained from Shuttle Radar Topography Mission SRTM website [20]. The DEM was processed and used to delineate the watershed and analyse the drainage patterns of the land surface terrain. Sub basin parameters such as slope gradient, slope length of the terrain, and the stream network characteristics such as channel slope, length, and width were derived from the DEM. Fig. 2 shows the DEM of the study area.

#### 2.3.2. Landuse and land cover

Land use is one of the most important factors that affect surface erosion, runoff, and evapo- transpiration in a watershed. Land use map of the study area was extracted from the Global Land Cover Characterization (GLCC) database and used to estimate vegetation and other parameters representing the watershed area. The GLCC database was developed by United States Geological Survey and has a spatial resolution of 1 km and 24 classes of land use representation [21].
Fig. 3 shows the land use and land cover of the study area while Table 1 provides information on land use of the study area. From Table 1, it could be inferred that the watershed is dominated by savannah type of vegetation having coverage of about 94% of the watershed.

#### 2.3.3. Digital soil data

SWAT requires different soil textural and physicochemical properties such as soil texture, available water content, hydraulic conductivity, bulk density and organic carbon content for different layers of each soil type. Digital soil data for the study was obtained from Harmonised Digital Soil Map of the world (HWSD v1.1) produced by Food and Agriculture Organization of the United Nations [22]. Table 2 provides detailed information on soil of the study area which shows that the area is dominated by clayey loam soil.

#### 2.3.4. Weather and flow data

SWAT requires daily meteorological data that can either be read from a measured data set or generated by a weather generator model. The weather data used in this study include daily precipitation, minimum and maximum temperature, relative humidity, wind, and solar radiation. These data were obtained from Nigerian Meteorological Agency (NIMET) Ilorin zone and covers a period of 10 y, i.e., 1 August 2004 to 31 July 2014. Observed flow data used for the calibration and validation of the model covers a period of 6 hydrological years i.e., August, 2008 to July, 2014.

### 2.4. Watershed Delineation and Creation of HRU

Watershed delineation tool on the Mapwindow GIS allows the user to sub-divide a watershed into several smaller hydrologically connected units which can be used for watershed characterization and modelling. The study area was divided into sub-basins using the GIS component of the SWAT model. The watershed was de lineated into 31 sub-basins which are then further subdivided into 42 HRUs that consist of homogeneous land use is to increase the accuracy of the prediction of model results from the sub basin level. Fig. 4 shows the stream network and the delineation of the watershed into sub-basins.

### 2.5. SWAT Setup and Configuration

The SWAT interface provides an option to choose the method to be adopted for the estimation of surface runoff (Curve Number or Green and Ampt method) and channel water routing (variable or Muskingum method). Also, for the estimation of potential evapo-transpiration (PET), methods such as Priestley, Penman-Monteith or Hargreaves are available. In this study, SWAT was run for a period of 10 hydrological years at a daily time step adopting Curve number method for estimating surface runoff from precipitation, variable storage method for channel water routing and Hargreaves method for estimating the PET from the watershed. The Hargreaves method was used for the estimation of PET because of its simplified equation which requires only two climatic parameters, that is, the temperature and incident radiation.

### 2.6. Performance Evaluation, Model Calibration and Validation

Performance of the SWAT model was assessed by comparing the models ability to match monthly values of observed flow (mean monthly discharge). In addition to comparing mean values for the calibration and validation periods, model performance was evaluated using statistical parameters such as the NSE and Coefficients of Determination (R2) as represented in Eq. (1) and (2), respectively.
##### (1)
$R2=[Σi(Qm,i-Qm¯)(Qs,j-Qs¯)]2Σm,j(Qm,j-Qm¯)2Σi(Qs,i-Qs¯)2$
##### (2)
$NSE=1-Σi(Qm,i-Qs,i)2Σi(Qm,i-Qm¯)2$
In Eq. (1) and (2), Qm is the measured discharge, Qs is the simulated discharge, $Qm¯$ is the average measured discharge and $Qs¯$ is the average simulated discharge. NSE values can range from 0 to 1. Perfect agreement between predicted and observed data results in NSE = 1; an NSE value of 0 indicates that the mean of the observed data is as accurate as the model predictions. For watershed scale modelling, NSE values of 0.36 to 0.75 are generally considered fair, while values greater than 0.75 indicate good model performance [23].

### 2.7. Computation of Theoretical Hydrokinetic Energy Potentials

Calibrated flow results of the model were used to estimate the theoretically available in-stream hydrokinetic energy potential of selected Rivers within the watershed. This was achieved using Eq. (3) in line with previous resource assessment of hydrokinetic energy potential [8].
##### (3)
$Pth=γ×Q×ΔH$
In Eq. (3), Pth is the theoretical hydrokinetic energy in Watts, Q is the segment specific average water discharge (m3/s), ΔH is the Hydraulic head (m) obtained from geo-processing of the DEM using the GIS component of the model and γ is the specific weight of water. The computations were finalised by exporting the data into Microsoft Excel spreadsheet software package to estimate the theoretical hydrokinetic energy potential along each of the River channels.
Fig. 5 shows the selected Rivers for the computation of hydrokinetic energy potential. Flow chart of the methodology adopted for the study is presented in Fig. 6.

### 3.1. Performance Evaluation, Model Calibration and Validation

Performance evaluation of the model was achieved by comparing the observed and simulated monthly inflow at Asa River for both the calibration and validation periods. In total, 14 parameters were selected to be calibrated through the Parasol optimization method. The model was calibrated with the observed monthly inflow from 2008 to 2010, and cross-validated with another set of independent data set from 2011–2014. Plot of observed data against calibrated results are presented in Fig. S1 and S2. Results showed a good correlation between the observed flow and the simulated flow, indicated by NSE and R2 of 0.76 and 0.85, respectively for calibration period, and NSE and R2 of 0.70 and 0.74, respectively for the validation period. Furthermore, most of the data for the calibration and validation period fall within the confidence interval of 95%. The plots of the standardised residuals for both the calibration and validation periods are as shown in Fig. S3 and S4, respectively.

### 3.2. Prediction of Average Daily and Annual Mean Discharge into Sub-basins

The results of average daily discharge into sub-basins of the watershed are as shown in Fig. S5. River reach in sub-basin 3 has the highest value of 2.79 m3/s while the reach in sub-basin 16 has the lowest value of 0.49 m3/s. Model results showed that sub-basin 3 has the highest value of 336.6 m3/s for maximum discharge in the watershed while sub-basin 16 has the lowest value of 38.91 m3/s for maximum discharge on the watershed. The results also showed that sub-basin 3 has the highest value of 1,053.4 m3/s for Annual Mean Discharge from the watershed. However, the lowest value of 16.02 m3/s was recorded in sub-basin 29.

### 3.3. Estimation of Hydrokinetic Energy Potential

In this study, theoretical hydrokinetic potential of Asa River, having a stream length of 76,713 m and an annual average flow of 59 m3/s was estimated as 41.63 MW. Oyun River has a length of 64,381m and an average annual stream flow of 125.5 m3/s was estimated as 98, 39 MW while River Awun, with an average annual flow of 1,053.4 m3/s and length of 31,736 m has the highest estimated hydrokinetic potential of about 154.82 MW. Details of theoretical hydrokinetic potentials of the Rivers can be found in Table S1. However, it is important to state that these values may be affected by factors such as turbine design, type, efficiency and mode of deployment and River use in the course of harnessing the potential on site.
The estimated hydrokinetic potentials most especially for Asa River in this work are comparable with the results obtained by [10] where the total theoretical hydrokinetic energy potential of selected basin upper catchment of Ikere Gorge Basin in Oyo State, Nigeria was analyzed for hydrokinetic energy and an estimated potential of 36.4 MW was obtained. Larger theoretical hydrokinetic potential estimated for River Awun (154.82 MW) may be due to the increased flow velocity of water at the location which emptied directly into River Niger, the largest River those transverses the country from Northern region to Southern region.

### 4. Conclusions and Future Outlook

Based on the outcome of this study, the following points can be concluded:
1. Naturally occurring hydrokinetic energy potentials exist along the three selected Rivers in the study area, namely, Asa, Oyun and Awun Rivers in Kwara State of Nigeria.

2. SWAT embedded in GIS environment can be employed in the assessment of hydrokinetic energy of rivers in other area of Nigeria and sub Saharan Africa.

3. Performance evaluation using coefficient of determination (R2) and NSE showed that the model performance is satisfactory in predicting the streamflow.

4. The estimate of the total annual theoretical hydrokinetic resource for Asa River was 42.63 MW, Oyun River was estimated as 98.39 MW while Awun River potential was obtained as 154.82 MW.

Therefore, with deployment of appropriate technology and commitment from the government and other stakeholders, estimated energy potentials along the selected Rivers could be harnessed in the nearest future and serve as alternate source of energy for the rural dwellers most especially settlements within the watershed that are not connected to national grid.

### Acknowledgment

Researchers acknowledge the support of the staff and students of the department of Civil Engineering, Kwara State University, Malete, Nigeria for their support towards the successful completion of this work.

### References

1. Bahleda M, Hosko MA. Assessment of water power potential and development needs. Electric Power Research Institute California Final Report. 2007.

2. Ladokun L, Ajao KR, Sule BF, Adeogun N. Hydrokinetic energy resource estimates of River ERO at Lafiagi, Kwara State, North-Central Nigeria. Nigerian J Technol Res. 2016;11:58–63.

3. Johnson JB, Pride DJ. River, tidal, and ocean current hydrokinetic energy technologies: Status and future opportunities in Alaska. Alaska energy authority by the Alaska Center for Energy and Power. 1November2010. Alaska:

4. Evan M. Hydrokinetic power: An analysis of its performance and potential in the Roza and Kittas canals [dissertation]. Olympia: Evergreen State College; 2012.

5. Mourad B, Enguerran G, Guillou S, Goossens X. Modelling river hydrokinetic energy in large rivers: Maroni river in French Guiana. 2016. cited 20 December 2017. Available from: https://hal.archives-ouvertes.fr/hal-01514240

6. Kalnacs A, Kalnacs J, Mutule A, Persis U. Methods for estimation of the river flow potential for hydrokinetic power generation. Latv J Phys Tech Sci. 2014;2:1–10.

7. Lalander E. Modelling hydrokinetic energy resource for in-stream energy converters [dissertation]. Uppsala: Uppsala Univ.; 2010. 1–57.

8. Canmet Energy. Assessment of Canada’s hydrokinetic power potential. Phase I Report methodology and data review. March2010. Ottawa:

9. Ladokun LL, Ajao KR, Adeogun AG. Hydrokinetic energy resource estimates of River ERO at Lafiagi, Kwara State, North-Central Nigeria. Nigerian J Technol Res. 11:58–63.

10. Wahab S, Ifabiyi IP, Adeogun AG. SWAT analysis of Ikere Gorge basin for hydrokinetic power estimation in selected rural settlements of Oke Ogun, Nigeria. Ruhuna J Sci. 2017;8:24–43.

11. Otuna JA, Onemano JI, Alayand AW. Assessment of the hydropower potential of Kangimi reservoir in Kaduna state Nigeria. Nigerian J Technol. 2012;31:300–307.

12. Vanguard Newspaper. NPC puts Nigeria’s population at 198 m [Internet]. C2018. [cited 12 December 2018]. Available from: https://www.vanguardngr.com/2018/04/npc-puts-nigerias-population-198m/

13. Emovon I, Samuel OD, Mgbemena CO, Adeyeri MK. Electric power generation crisis in Nigeria: A review of causes and solutions. Int J Integr Eng. 2018;10:47–56.

14. Federal Ministry of Water Resource. Dams and reservoir operations [Internet]. [cited 2 December 2017]. Available from: http://waterresources.gov.ng/dams-and-reservoir-operations/

15. Ladokun LL, Ajao KR, Sule BF. Hydrokinetic energy conversion systems: Prospects and challenges in Nigerian hydrological setting. Nigerian J Technol. 2013;32:371–378.

16. Jimoh HI, Ishola AL. Problems of suspended sediment loads in Asa River catchment, Ilorin, Nigeria. Pakistan J Soc Sci. 2009;6:19–25.

17. Arnold JG, Williams JR, Maidment DR. Continuous-time water and sediment-routing model for large basins. J Hydraul Eng. 1995;121:171–183.

18. Neitsch SL, Arnold JG, Kniry JR, Scrinivasan R, Williams JR. Soil and water assessment tool–user’s manual–version 2000. Temple (TX): Blackland Research Centre; 2002.

19. Neitsch SL, Arnold JG, Kiniry JR, Williams JR. Soil and water assessment tool–theoretical documentation–version 2005. Temple (TX): Blackland Research Centre; 2005.

20. CGIAR. SRTM 90m digital elevation data [Internet]. [cited 20 July 2017]. Available from: http://srtm.csi.cgiar.org/

21. USGS. Global land cover characterization [Internet]. [cited 25 November 2017]. Available from: https://www.usgs.gov/centers/eros/science/usgs-eros-archive-land-cover-products-global-land-cover-characterization-glcc?qt-science_center_objects=0#qt-science_center_objects

22. Nachtergaele F, Velthuizen HV, Verelst L. Harmonized world soil database. In : Fod and Agriculture Organization of the United Nations; March 2009; Rome, Italy. 2008. Laxenburg: IIASA; 2008

23. Motovilov YG, Gottschalk L, England K, Rodhe A. Validation of distributed hydrological model against spatial observations. Agric Forest Meteorol. 1999;98–99:257–277.

##### Fig. 1
Map of Nigeria showing the location of the watershed.
##### Fig. 2
Digital elevation model of the study area attributed with stream network.
##### Fig. 3
Landuse map of the study area.
##### Fig. 4
Delineation of watershed into 31 sub-basins and 42 hydrological response units (HRUs).
##### Fig. 5
Selected Rivers for the computation of hydrokinetic energy potential.
##### Fig. 6
Flow chart of the stages involved in the estimation of hydrokinetic energy of the selected rivers.
##### Table 1
Information on Land use of the Study Area
S/N SWAT Code Description Area(Ha) % of Watershed
1 SAVA Savannah 603,710.23 94.44
2 CRWO Cropland/woodland mosaic 10,288.11 1.61
3 CRDY Dryland cropland and pasture 56.39 0.01
4 SHRB Shrubland 749.11 0.12
5 CRGR Cropland/grassland mosaic 749.11 0.12
6 WATB Waterbodies 3,182.47 0.50
7 BSVG Barren or sparsely vegetated 589.19 0.09
##### Table 2
Information on Soil Texture of the Study Area
S/N SWAT Code Description Area(Ha) % of Watershed
1 I-Ql-Re-1281 Loam 32,824.66 5.13
2 Lf1-1421 Sandy loam 5,526.61 0.86
3 Lf12-a-1430 Clayey loam 159,378.61 24.93
4 Lf26-a-1443 Sandy clayey loam 415,467.82 64.99
5 Lf49-1475 Sandy clayey loam 6,126.74 0.96
6 Nd3-1565 Sandy loam 19,914.63 3.13
Total 639,239.07 100.00
TOOLS
Full text via DOI
Supplement
E-Mail
Print
Share:
METRICS
 1 Crossref
 1 Scopus
 1,866 View


text:
# General strategies to write big code in Mathematica?

I think after six months of exposure to Mathematica and the Wolfram Language I am fairly OK with writing short codes and snippets. However, what are some general strategies to use in order to write big blocks of code?

• If you can decompose your "big block" into a bunch of subroutines, it should help a great deal with taming the complexity. – J. M. will be back soon Mar 13 '16 at 7:37
• I don't agree with the people who are trying to close this as a duplicate. If we don't have a question that specifically addresses this question, I say let it be. – C. E. Mar 13 '16 at 11:12
• @C.E. Same here. To the closers is there any good reason to close? I personally don't think so. This is an important topic, which, AFAIK, wasn't discussed here as a separate one yet. The linked question does not answer this. – Leonid Shifrin Mar 13 '16 at 11:21
• @JeelShah I think Wolfram Language is a compact and elegant language. Last year, I use Mathematica to design a CAGD package with the help of the classical book "The NURBS Book". The author used approximately $300,000+$ lines C code to implement a full NURBS package. However, the Mathematica code that written in my package is only $5,000+$ lines. Lastly, owing to Mathematica is a interpreter language, rather than a compiled language, so generally the speed/performance of functions that written inMathematica language is slower than the function written in C. – xyz Mar 14 '16 at 2:26
• @ShutaoTANG Actually the fun starts when you write Mathematica code in such a way, using Compile (to C) and LibraryLink, that you gain back a lot of speed. I like that in this way I can get quite nice generated C-code, for, say, nested loops, but keep the rest in high level Mathematica with all its nice cosmetics. – Rolf Mertig Mar 14 '16 at 21:10

# Managing the complexity, II: controlling complexity on the smaller scale

There are a few things you can do to control and reduce the complexity of your code, even on the small scale - long before you move to packages and split code into several files.

## Effective use of the core data structures

This is probably the first thing to mention. The most important core data structures are Lists and Associations. Mastering them and using them effectively goes a long way towards writing much better Mathematica code.

Some of the properties which make both Lists and Associations so effective are:

• They are very well integrated into the language
• They are polymorphic data structures. Lists can hold elements of any type, and Associations can use elements of any type both for keys and for values.
• They are very universal. In particular, Lists can be used for arrays, sets, trees, etc., and Associations implement a very general key - value mapping abstraction.
• They are mostly immutable, with a very limited form of mutability. This is, in fact, a big advantage.
• Using them with a functional programming style leads to very compact code doing non-trivial data transformations fast. This both reduces the code bloat and improves the code speed.
• They offer a fast and cheap way to do exploratory programming and prototyping, where you don't have to create new types, so can create and change complex data structures on the fly.

However, in the long term, one has to be aware of certain flaws as well:

• Lists:
• Adding and removing elements is O(n) operation, where n is the length of the list
• Associations
• Are rather memory-hungry
• Element by element modifications can be rather slow. Even though Associations themselves have roughly O(1) complexity for these operations, one still have to do a top-level iteration to, for example, build an association element by element, and top-level iteration is slow. In other words, there is no analogue of packed arrays for associations. In some cases, one can use functions like AssociationThread, which operate on many keys and values at once.
• Common
• It is easy to get regression bugs from changes in code, due to weak typing

## Code granularity

In most cases, it is much better to split your code in a number of fairly small functions, each one doing some very specific task. Here are a few suggestions regarding that:

• Use small functions (just a few lines of code each)
• Avoid side effects as much as possible
• In particular, prefer With to Module when possible
• Write code in a style that promotes function composition
• Use operator forms and currying (available for a number of built-in functions since V10)

### Example: simplistic DOM viewer

Below is the code of a rudimentary viewer for a DOM structure of an HTML page:

ClearAll[view,processElement,color,grid, removeWhitespace, framed];

framed[col_] := Framed[#, Background -> col] &;

color[info_] :=
Replace[
info,
s_String :>
Mouseover[framed[LightYellow][s], framed[LightGreen][s]],
{1}
];

removeWhitespace[info_]:=
DeleteCases[info,s_String/;StringMatchQ[s,Whitespace]];

grid[info_List]:=Grid[{info},ItemSize->Full,Frame-> All];

processElement[tag_,info_]:=
OpenerView[{tag, grid @ color @ removeWhitespace @ info}];

view[path_String]:=
With[{expression = Import[path,"XMLObject"]},
view[expression] /; expression =!=$Failed ]; view[expr_]:= Framed[ expr//.XMLObject[_][_,info_,_]:>info//. XMLElement[tag_,_,info_]:>processElement[tag,info] ];  You can call view[url-from-the-web] to view some web page,for example view["https://www.wikipedia.org/"]  This is what I call granular code: it contains a few really tiny functions, which are very easy to understand and debug. ### Example: modeling and visualizing random walks This one was a real question asked by someone in the Russian-speaking Mathematica online group. It is valuable since this is a real problem, and it was originally formulated in a procedural style. The problem is to model a 2-dimensional random walk with certain step probabilities, which are constants (don't depend on the previous steps). The question asked is to find a probability to return to the point of origin in less than a given number of steps. This is done using essentially Monte-Carlo simulation, running the single walk simulation multiple times, and finding how many steps it took to return back, for a particular experiment. Here is the original code. The problem settings (I keep the original code): (* v - the possible point's displacements *) v = {{1, 1}, {-1, 1}, {-1, -1}, {1, -1}}; (* p - probabilities for the next step *) p = Array[{(#)/4, #} &, {4}]; (* A number of experiments *) n = 23; (* Maximal number of steps *) q = 100; (* Number of repetitions for all experiments, used for computation of a mean and standard deviation for the probability density for return *) Z = 500;  The actual computation (* The choice of the next step *) Step[R_] := v[[Select[p, #[[1]] >= R &, 1][[1]][[2]]]]; (* Array initialization. m[[i]] gives a number of successful returns in i-th run *) m = Array[#*0 &, {Z}]; (* Running the experiments *) For[k=0,k<Z,k++ For[j=0,j<n,j++, (* Initial position of a point *) X0={0,0}; (* Making the first step *) i=1; X0+=Step[RandomReal[]]; (* Move until we return to the origin, or run out of steps *) While[(X0!={0,0})&&(i<q),{X0+=Step[RandomReal[]],i++}]; (* If the point returned to the origin, increment success counter *) If[X0=={0,0},m[[k]]++]; ]; ];//AbsoluteTiming (* {5.336, Null} *)  Here is the visualization of the experiment (basically, the unnormalized empirical CDF and PDF): GraphicsGrid[{{ ListPlot[Sort[m/n], Joined -> True, PlotRange -> All, PlotStyle -> Thick, Filling -> Axis], Histogram[m/n, Automatic] }}]  The original complaint was that the code was slow, and the asker was looking for ways to speed it up. Here is my version of the code, rewritten in a functional granular style: ClearAll[returnedQ,randomSteps,runExperiment,allReturns,stats]; (* Determine if the point returned*) returnedQ[v_,steps_]:=MemberQ[Accumulate[v[[steps]]],{0,0}]; (* Generate random steps *) randomSteps[probs_,q_]:=RandomChoice[probs->Range[Length[probs]],q]; (* Run a single experiment *) runExperiment[v_,probs_,q_]:= returnedQ[v,randomSteps[probs,q]]; (* Run all n experiments *) allReturns[n_,q_,v_,probs_]:= Total @ Boole @ Table[runExperiment[v,probs,q],{n}] (* Gather all statistics from z runs *) stats[z_,n_,q_,v_,probs_]:=Table[allReturns[n,q,v,probs],{z}];  We run it as (m = stats[Z, n, q, v, {1/4, 1/4, 1/4, 1/4}]); // AbsoluteTiming (* {0.411055, Null} *)  It ends up 10 times faster, but also the above code is, at least for me, much more readable - and you can easily test all individual functions, since they don't depend on anything that has not been passed to them explicitly. Now, here is a version of the same code, expressed as a one-liner: statsOneLiner[z_,n_,q_,v_,probs_]:= Table[ Total @ Boole @ Table[ MemberQ[ Accumulate[v[[RandomChoice[probs->Range[Length[probs]],q]]]], {0,0} ], {n} ], {z} ];  What I would say is that I strongly prefer the granular version, in all cases but those where the condensed one offers far superior performance, and only if this is critical for the problem. In this particular case, the performance is the same, and in most other cases it also won't be worth it to keep such code, since it is much harder to understand. In any case, to me this example serves as another good illustration of the advantages and superiority of functional programming done in a granular fashion, and I hope it additionally illustrates my point about the importance of granularity. ## Function composition Writing code in this style is very beneficial for readability, extensibility and the ease of debugging. Do it, when you can. ### Example: inverting many to many relationships I will borrow this one from this answer. The function below inverts many-to-many relationship encoded in an association: ClearAll[invertManyToMany]; invertManyToMany[start_Association]:= Composition[ Merge[Identity], Reverse[#,{2}]&, Catenate, Map[Thread], Normal ] @ start;  Here is an example of use: invertManyToMany @ Association[{ "programming" -> {1, 2, 3},"graphics" -> {2, 3, 4}, "plotting" -> {1, 3, 4, 5}} ] (* <| 1 -> {"programming", "plotting"}, 2 -> {"programming", "graphics"}, 3 -> {"programming", "graphics", "plotting"}, 4 -> {"graphics", "plotting"}, 5 -> {"plotting"} |> *)  But here I just want to stress the way the function is written: using Composition and operator forms makes the code much more transparent and much easier to debug and extend. To debug, you basically need to stick something like showIt in between any two transformations in the chain, and to extend, you can simply add transformations. ## Function overloading When you define functions using patterns, you can use function overloading - giving several definitions to a single function, on various number / types of arguments. Languages which support overloading, have mechanisms for automatic dispatch to the right definition, given specific input arguments. This automation can be used to simplify programmer's life and write more expressive code. Mathematica fully supports overloading via its core pattern-matching engine, and in fact its pattern-matching capabilities can be thought of as "overloading on steroids" in this context, compared to other languages - even those which support multiple dispatch. You can actually often design your code in a such a way as to maximally utilize this option. Functions written in such a style are typically (not always though) much more readable and extensible, than if you would have a single large Switch (or worse yet, nested If) inside the body. The reason has partly to do with the fact that this technique is roughly equivalent to introduction of ad-hoc mini type systems (since you overload on function arguments, which you check using patterns, and thus defining weak types for them), and partly because Mathematica allows multiple dispatch, which is much more powerful than a single - argument dispatch available in many other languages. I will illustrate this with a single example taken from the RLink module source code: this single function determines the type of all RLink objects, either sent to R from Mathematica, or received from R: ClearAll[typeOf]; typeOf[o_?outTypeInstanceQ]:=o@getType[]@getStringType[]; typeOf[var_String]:=typeOf[var,getRExecutor[]]; typeOf[var_String, RExecutor[exec_]]:= With[{type=exec@getRObjectType[var]}, type/;type=!=Null ]; typeOf[var_String,_RExecutor]:="UnknownRCodeType"; typeOf[RVector[type_,data_,att_]]:=type; typeOf[RList[data_,att_]]:="list"; typeOf[RNull[]]:="NULL"; typeOf[_]:=Throw[$Failed,error[typeOf]];

This example illustrates two more quite useful tricks: use local variables shared between the body and the condition of the rule, and use the catch-all pattern to throw local (internal) exception - but these I will discuss separately.

To summarize advantages of this method:

• Easier to read, understand, write, and debug such code
• Code written in this way is more extensible
• Often you can get rid of intermediate variables, that would be necessary otherwise

Some of the things to watch for are, though:

• You have to keep an eye on the relative generality of definitions
• In some rare cases, you may need to manually reorder definitions
• You can't use Compile on functions which use rules and patterns

## Small scale encapsulation: inner functions

This is a form of encapsulation, where you introduce inner functions, local to the Module, Block, or With scoping constructs that you use to encapsulate your local variables / state. The advantage of this technique is that you can achieve a better level of modularity and readability of your code on a smaller scale, without using such a heavy tool as contexts and packages.

### Example: directory traversals with skips

Here is an example I took from this Mathgroup post and modified:

ClearAll[clearSkip,setSkip,dtraverse, shallowTraverse];

shallowTraverse[dir_String,dirF_,fileF_]:=
Scan[
If[FileType[#]===Directory,dirF[#],fileF[#]]&,
FileNames["*",dir]
];

Module[{skip},
clearSkip[]:=(Clear[skip];skip[_]=False);
setSkip[dir_String]:=skip[dir]=True;
dtraverse[dir_String,dirF_,fileF_]:=
Module[{traverse,level=0, withLevel},
withLevel = Function[fn, level++;fn[level];level--];
traverse[ldir_]:=
withLevel[
Function[lev,
dirF[ldir,lev];
If[!TrueQ[skip[ldir]],
shallowTraverse[ldir,traverse, fileF]
];
]];
shallowTraverse[dir,traverse,fileF]
]
];

The usage is:

dtraverse[directory,dirF,fileF],

where dirF accepts 2 parameters: subdirName,level, and fileF accepts 1 parameter - file name. You can use this to traverse a directory tree, applying arbitrary functions to files and directories at a specified level, and can set at runtime directories which have to be skipped entirely.

Before we run this code, a few words about it. It is all built on inner functions and closures. Note that all of the clearSkip, setSkip and dtraverse are closed over a local variable skip. Moreover, withLevel and traverse are inner closures, closed over level and skip, fileF and dirF, respectively. What do I buy with closures? Better composition, and better code structuring. Because I don't have to explicitly pass the parameters, I can, for example, pass traverse directly as a parameter to shallowTraverse, making the code easier to read and understand.

The code structure here is very transparent. I view nested directory traversal with functions fileF and dirF as a shallow traversal, where fileF gets applied to files, while to the sub-directories we apply the traverse function. Now, what do I buy with factoring out withLevel? I could've easily wrapped level++;code;level-- in the body of traverse. The answer is, I separate the side effect. Now I could test the inner Function[lev, ...] in isolation, at least in principle.

Let us now see what the run-time skip facility can give us. Here I will run through the entire directory tree for the $InstallationDirectory, but only collect the names of the first-level subdirectories: clearSkip[]; Map[FileBaseName] @ Reap[ dtraverse[$InstallationDirectory, If[#2 === 1, Sow[#1]] &, 1 &]
][[2, 1]] // AbsoluteTiming

(*
"Frameworks", "Library", "MacOS", "Resources", "SystemFiles"}}
*)

Now I do the same, but instruct the code to skip the inner sub-directories' trees, setting setSkip at runtime:

clearSkip[];
Map[FileBaseName] @ Reap[
dtraverse[\$InstallationDirectory,If[#2 === 1, Sow[#1]; setSkip[#1]] &, 1 &]
][[2, 1]] // AbsoluteTiming

(*
"Frameworks", "Library", "MacOS", "Resources", "SystemFiles"}}
*)

And get the same result, only 1000 times faster. I think this is pretty cool given that it only took a dozen lines of code to implement that. And using inner functions and closures made the code clear and modular even on such a small scale, and allowed to cleanly separate state and behavior.

### Example: Peter Norvig's spelling corrector in Mathematica

In the following example this idea is pushed to the extreme. Here is where it comes from. It is hard to beat the clarity and expressiveness of Python, but at least I gave it a try.

Here is a training data (it takes some time to load this):

text = Import["http://norvig.com/big.txt", "Text"];

Here is the code I ended up with (I cheated a bit by abbreviating a number of built-ins using With, because in the original post of Norvig, there was a kind of competition between languages, and I wanted the code as short as possible, without losing readability. But I ended up liking it):

With[{(* short command proxies *)
fn=Function,join=StringJoin,lower=ToLowerCase,rev=Reverse,
smatches=StringCases,wchar=LetterCharacter,chars=Characters,
inter = Intersection,dv = DownValues,len=Length,
flat=Flatten,clr=Clear,replace=ReplacePart,hp=HoldPattern
},
(* body *)
Module[{rlen,listOr,keys,words,model,train,clrmemo,
transp,del,max,ins,repl,edits1,known,knownEdits2
},
(* some helper functions - compensating for a lack of built-ins *)
rlen[x_List]:=Range[len[x]];
listOr=fn[Null,Scan[If[#=!={},Return[#]]&,Hold[##]],HoldAll];
keys[hash_]:=keys[hash]=Union[Most[dv[hash][[All,1,1,1]]]];
clrmemo[hash_]:= If[UnsameQ[#/.Most[dv[keys]],#]&@hp[keys[hash]],keys[hash]=.];
(* Main functionality *)
words[text_String]:=lower[smatches[text,wchar..]];
With[{m=model},
train[features_List]:=
(clr[m];clrmemo[m];m[_]=1;Map[m[#]++&,features];m)
];
transp[x_List]:=
Table[replace[x,x,{{i},{i+1}},{{i+1},{i}}],{i,len[x]-1}];
del[x_List]:=Table[Delete[x,i],{i,len[x]}];
retrain[text_]:=
With[{nwords=train[words[text]],alphabet =CharacterRange["a","z"]},
ins[x_List]:=flat[Outer[Insert[x,##]&,alphabet,rlen[x]+1],1];
repl[x_List]:=flat[Outer[replace[x,##]&,alphabet,rlen[x]],1];
max[set:{__String}]:=Sort[Map[{nwords[#],#}&,set]][[-1,-1]];
known[words_]:=inter[words,keys[nwords]]
];
Attributes[edits1]={Listable};
edits1[word_String]:=
join @@@ flat[Through[{transp,ins,del,repl}[chars[word]]],1];
knownEdits2[word_String]:=known[flat[Nest[edits1,word,2]]];
(* Initialization *)
retrain[text];
(* The main public function *)
correct[word_String]:=
max[listOr[known[{word}],known[edits1[word]],knownEdits2[word],{word}]];
]
];

Here are some tests:

correct["proed"] // Timing

(* {0.115998, "proved"} *)

correct["luoky"] // Timing

(* {0.01387, "lucky"} *)

The above code illustrates another thing about inner functions: you may use them also to significantly change the way the code looks inside Module. Of course, we could write it all in a single incomprehensible one-liner, but that would hardly benefit anyone.

### Summary

I personally use inner functions all the time, and consider them an important tool for improving small-scale encapsulation, structure, readability and robustness of the code.

One thing to watch out for is that in some cases, inner functions are not garbage - collected automatically. This usually happens when some external objects point to them at the time when they are defined. This may or may not be acceptable, depending on your circumstances. There are also ways to avoid it, such as using pure functions (which, however, can't be easily overloaded and are generally less expressive since you can't easily do pattern-based arguments destructuring and tests for them).

## Final example: Huffman encoding

To illustrate many of the points I mentioned above, I will here provide my re-implementation of the Huffman encoding algorithm, based on the code from David Wagner's excellent book. So I refer to his exposition for details on the algorithm and ideas used. I rewrote it to use Associations, and made it purely functional, so that there is no mutable state involved whatsoever, anywhere in the code.

msg = "she sells sea shells by the sea shore";
chars = Characters@msg

(*
{"s", "h", "e", " ", "s", "e", "l", "l", "s", " ", "s", "e", "a",
" ", "s", "h", "e", "l", "l", "s", " ", "b", "y", " ", "t", "h",
"e", " ", "s", "e", "a", " ", "s", "h", "o", "r", "e"}
*)

### Building Huffman tree

Here is all the code needed to build Huffman tree from an arbitrary list of elements:

ClearAll[combine];
combine[{x_}]:={x};
combine[{{n_,x_},{m_,y_},rest___}]:= Union[{{m+n,{x,y}},rest}]

huffmanTree[elems_List]:=
With[{pool = Sort @ Reverse[#,{2}]& @ Tally @ elems},
FixedPoint[combine,pool][[1,2]]
];

Here is the tree in our case:

tree = huffmanTree@chars

(*
{{{{"y", "a"}, "h"}, "s"}, {{"l", {{"b", "o"}, {"r", "t"}}}, {" ", "e"}}}
*)

### Encoding

Here is all the code needed to encode the message, given a Huffman tree:

ClearAll[huffEncode];
huffEncode[chars_List, tree_]:=
Composition[
Flatten,
Lookup[chars],
AssociationMap[First[Position[tree,#]-1]&],
Union
] @ chars

Now we encode our message:

encoded  = huffEncode[chars, tree]

(*
{0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0,
0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1,
1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0,
0, 1, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1,
1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1,
0, 1, 1, 1}
*)

### Decoding

Here is the code to decode the message:

ClearAll[extract];
extract[tree_][{subtree_, _}, code_Integer]:= extract[tree][subtree[[code+1]]];
extract[tree_][elem_String]:={tree, elem};
extract[tree_][subtree_]:={subtree, Null};

ClearAll[huffDecode];
huffDecode[tree_, codes_]:=
DeleteCases[Null] @ FoldList[extract[tree],{tree,Null},codes][[All,2]];

So that we have

StringJoin @ huffDecode[tree, encoded]

(* "she sells sea shells by the sea shore" *)

### Notes

I think, this example illustrates very well the kind of economy and simplicity that is possible to get from a combination of functional programming, very granular code, function overloading, function composition, operator forms / currying (note that I actually introduced currying / operator form also for the user-defined extract function), and the user of core data structures (Lists and Associations) in Mathematica.

All code contains absolutely no mutable state. Except for inner functions, it uses all of the techniques I described above. The result is a tiny program that solves a non-trivial problem, and while there is no room here for code dissection, it is very easy to take this code apart and understand what goes on at each step. In fact, it is mostly clear how it works just from looking at the code.

Of course, the main credit goes to David Wagner, I just made a few changes to utilize a few recent additions like Associations and completely remove any mutable state.

• This shouldn't be a community wiki... I think you deserve all the points that come with all these answers! Nobody else is going to be able to add to this answer authoritatively, so it really isn't a "wiki"... :D If you would agree, I'll revert the change :) – rm -rf Mar 14 '16 at 15:33
• @R.M. Well, I didn't want to multiply the points I get for the sub-answers, but OTOH it might be good to retain the ownership indeed, since I want to preserve the posts mostly in the way I have in mind (I never object to corrections and suggestions, so will anyway incorporate those which make sense). So well, I guess I'd agree with your suggestion. Thanks :) – Leonid Shifrin Mar 14 '16 at 15:38
• Heads up, I am taking a counterpoint to one of your examples, and as always I look forward to your rebuttal: (110067) – Mr.Wizard Mar 15 '16 at 10:31
• It also seems that v[[RandomChoice[probs -> Range[Length[probs]], nn]]] can be simplified to RandomChoice[probs -> v, nn] . If there is a problem with packed arrays, it can be resolved by making v a packed array. Perhaps this solution did not work in previous versions. (@Mr.Wizard) – Jacob Akkerboom Mar 15 '16 at 18:39
• @JacobAkkerboom I am most happy to provide them. This is an important topic, and I am surprised that this question hasn't been asked long time ago. I do hope to find time to work some more on the answers, since there still is a bunch of stuff I consider useful which I did not add here yet. – Leonid Shifrin Mar 15 '16 at 20:40

# Preamble

I had a talk devoted specifically to this topic, on Second Russian WTC in 2014. Unfortunately, it is in Russian. But I will try to summarize it here.

Since this post is becoming too long, I decided to split it to several smaller ones, each dedicated to some particular set of methods / techniques. This one will contain a general / conceptual overview. As I add more specific parts, the links to them will be added right below this line.

• Controlling complexity on the smaller scale
• Effective use of core data structures
• Code granularity
• Function composition
• Small-scale encapsulation, scoping, inner functions
• Using powerful abstractions
• Higher - order functions
• Closures
• Abstract data types, stronger typing
• Macros and other metaprogramming techniques (to be added, not yet there)

### Problems

From the bird's eye perspective, here is a list of problems typically associated with the large-scale development, which is largely language-agnostic:

• Strong coupling between the modules
• Losing control over the code base as it grows (it gets too hard to keep in mind all things at once)
• Debugging becomes harder as the code base grows
• Loss of flexibility for the project, it becomes harder to evolve and extend it

### Techniques

Some of the well-known methods to tame projects complexity include:

• Separation between interfaces and implementations, ADTs
• Module system / packages / namespaces
• Mechanisms of encapsulation and scoping, controlling visibility of functions and variables
• Stronger typing with type inference
• Use of more powerful abstractions (if the language allows that)
• Object models (in object - oriented languages)
• Design patterns (particularly popular in OO - languages)
• Metaprogramming, code generation, DSLs

### Goals

All these methods basically help to reach a single goal: improve the modularity of the code. Modularity is the only way to reduce complexity.

To improve modularity, one usually tries to:

• Separate the general from the specific

This is one of (if not the) most important principles. It is always good to separate degrees of abstraction, since mixing them all together is one of the major sources of complexity.

• Split code into parts (functions, modules, etc)

Note that splitting into parts doesn't mean a trivial splitting of code into several parts - this helps only a little. It is a more complex process, where one has to identify a set of general abstractions, a mix of which, parametrized with the specifics, results in the most economical and simple implementation - and then separate those abstractions into separate functions or modules, so that the remaining specific part of the code is as small and simple as possible.

This may be non-trivial, because it may not be apparent from the specific implementation, what these parts are, and one has to first find points of generalization, where the code must first be generalized - then those "joints" will become visible. It takes some practice and a habit of thinking this way, to easily identify these points.

It also greatly depends on what's in your tool set: these "split points" will be different for say, procedural, functional and OO programming styles, and that will result in different splitting patterns and at the end, different code with different degrees of modularity.

• Decrease coupling between the parts

This basically means decreasing their inter-dependencies, and have well-defined and simple interfaces for their interaction. To reach this goal, one frequently adds levels of indirection, and late decision - making.

• Increase the cohesion for each part (so that the part is much more than the sum of its sub-parts)

This basically means that the components of each part don't make much sense taken separately, much like parts of the car's engine - you can't take out any one of them, they are all inter-related and necessary.

• Make decisions as late as possible

A good example in the context of Mathematica is using Apply: this postpones the decision on which function is called with a given set of arguments, from write-time to run-time.

Late decision - making decreases coupling, because interacting parts of the code need less information about each other ahead of time, and more of that information is supplied at run-time.

# General things

Here I will list some general techniques, which are largely language - agnostic, but which work perfectly well in Mathematica.

## Embrace functional programming and immutability

A lot of problems with large code bases happen when the code is written is stateful style, and state gets mixed with behavior. This makes it hard to test and debug separate parts of the code in isolation, since they become dependent on the global state of the system.

Functional programming offers an alternative: program evaluation becomes a series of function applications, where functions transform immutable data structures. The difference in resulting code complexity becomes qualitative and truly dramatic, when this principle is followed down to the smallest pieces of code. The key reason for this is that purely functional code is much more composable, and thus much easier to take apart, change and evolve. To quote John Hughes ("Why functional programming matters"),

"The ways in which one can divide up the original problem depend directly on the ways in which one can glue solutions together. "

I actually highly recommend to read the entire article.

In Mathematica, the preferred programming paradigms, for which the language is optimized, are rule-based and functional. So, the sooner one stops using imperative procedural programming and moves to functional programming in Mathematica, the better off one will be.

## Separate interfaces and implementations

This has many faces. Using package and contexts is just one, and rather heavy, way to do that. There exist also ways to do that on the smaller scale, such as

• Creating stronger types
• Using the so-called i - functions
• Inserting pre and post-conditions in functions

## Master scoping constructs and enforce encapsulation

Mastering scoping is essential for scaling to larger code bases. Scoping provide a mechanism for information - hiding and encapsulation. This is essential for reducing the complexity of the code. In non-trivial cases, it is quite often that, to achieve the right code structure, even inside a single function, one may need three, four or even more levels of nesting of various scoping constructs (Module, Block, With, Function, RuleDelayed) - and to do that correctly, one has to know exactly what the rules of their mutual interaction are, and how to bend those rules if necessary. I can't overemphasize the importance of scoping in this context.

## Separate orthogonal components in your code

This is a very important technique. It often requires certain advanced abstractions, such as higher-order functions and closures. Also, it requires some experience and certain way of thinking, because frequently code doesn't look like it can be factored - because for that certain parts of it should be rewritten in a more general way, yet it can be done. I will give one example of this below, in the section on higher-order functions.

## Use powerful abstractions

Here I will list a few which are particularly useful

• Higher-order functions
• Closures
• Function composition
• Strong types
• Macros and other meta-programming devices

## Use effective error-reporting in internal code, make your code self-debugging

There are a number of ways to achieve that, such as

• Using Assert
• Setting pre and post - conditions
• Throwing internal exceptions

All them combined, lead to a much simpler error diagnostics and debugging, and also greatly reduce regression bugs

## Use unit tests

There has been enough said about the usefulness of unit tests. I just want to stress a few additional things.

• Mathematica meta-programming capabilities make it possible and relatively easy to simplify generation of such tests.
• The extremely fast development cycle for the prototyping stage somewhat flies in the face of unit-testing, since code changes so fast that writing unit tests becomes a burden. I would recommend to write them once you move from a prototype to a more stable version of a particular part of your code.

# Topics not covered yet (work in progress)

To avoid making this post completely unreadable, I did not cover a number of topics which logically belong here. Here is an incomplete list of those:

• More details about packages and contexts
• Error reporting and debugging
• Using metaprogramming, macros and dynamic environments
• Using development tools: Workbench, version control systems
• Some advanced tools like parametrized interfaces

# Summary

There are a number of techniques which may be used to improve the control over code bases as they grow larger. I tried to list a few of them and give some examples to illustrate their utility. These techniques can be roughly divided into a few (overlapping) groups:

• Small-scale techniques
• Effective use of core data structures
• Code granularity
• Small-scale encapsulation, scoping, inner functions
• Function composition
• Large-scale techniques
• Packages and contexts
• Factoring orthogonal components
• Separation of interfaces and implementations
• Using powerful abstractions
• Abstract data types, stronger typing
• Closures
• Higher - order functions
• Macros and other metaprogramming techniques

This is surely not an ideal classification. I will try to make this post a work in progress and refine it in the future. Comments and suggestions more than welcome!

• Did you say "summarise"? :) – Martin Ender Mar 13 '16 at 17:17
• @MartinBüttner This is true. That was a 1.5 hour talk with about 10x more stuff in there than is here in this post. But I will probably remove some of the examples or replace them with smaller ones. It would take another several hours to make it qualitatively better, and I don't have more time for this today. Consider this an alpha version of the post :) – Leonid Shifrin Mar 13 '16 at 17:21
• Don't worry, I was kidding. It's a great answer. :) – Martin Ender Mar 13 '16 at 17:22
• @LeonidShifrin I watched your linked presentation. My impression is that in 2014 you were studying Japanese and you switched the language localization of your OS to Japanese in order to practice reading Kana. – Anton Antonov Mar 13 '16 at 17:54
• @bobthechemist Using Module in most cases means that you intend to modify the varables in the body of it, so with Module you introduce local mutable state. Using With means that you can't modify your local variables (actually, constants) in its body, so your code is completely side effect - free (unless there are some additional mutable variables in your code) , which is generally better. This means, in particular, that you can extract any part of the code inside With and test it separately - which you can't generally do with Module. – Leonid Shifrin Mar 14 '16 at 7:32

Here are some advices from my experience.

• Explore new ideas with the Mathematica frontend. Don't hesitate to use sections and subsections in the frontend to structure your work and experiment various possibilities.
• When you have instructions that work, package them into functions, still in the frontend.

It's practical to select all the useful instruction cells while holding Ctrl (in Windows), to copy them somewhere in the notebook so that they are following each other and to merge them. Then you just have to add a Module around the code, localize variables and add a function declaration with arguments.

• Then package them into packages. I do it like explained here. You can also do it from a notebook.
• Use Wolfram Workbench. It's really important from my point of view for big projects as having a debugger is very important. Also you can rename variables across multiple packages (files) which is very convenient.

(The only non obvious thing in using Workbench, is when you're at a breakpoint in debug mode and want to abort the evaluation, you need to abort the evaluation in the Mathematica frontend using Ctrl+. for example and continue the evaluation in Workbench.)

• Once you already have a project big enough, you can write some functions directly in Workbench.
• Write unit tests, before or right after writing a new code that works. Workbench handles unit tests.
• Use code versioning, for example Git with the plugin Egit in Eclipse (that you will use if you use Wolfram Workbench).
• Reuse, reuse, reuse. Never write twice the same thing.
• Use the function Echo or this more sophisticated utility to print values from deep inside your code.

# Managing the complexity III: using powerful abstractions

In this section I will list a few techniques which allow one to write more modular code and better separate the concerns, by using certain powerful abstractions provided by or possible to have in the Mathematica.

## Higher-order functions

These are functions which take other functions as arguments. In Mathematica, a number of core built-in functions like Map and Apply are higher-order functions.

The utility of this construct can be seen most clearly within the functional programming paradigm. Higher-order functions can be used to parametrize generic functionality, where custom behavior is injected with functional arguments. This allows one to easily separate generic functionality from the specific.

### Trivial example: Select

One trivial example of a built-in higher-order function is Select. We can write a more specific version of Select that would select numbers larger than a threshold by parametrizing Select with an appropriate test function:

ClearAll[mySelect];
mySelect[l_List, threshold_]:=Select[l,#>threshold&]

So that

mySelect[Range[10], 5]

(* {6, 7, 8, 9, 10} *)

Note that the test function #>threshold& is in fact a closure, closed over threshold and created at run-time.

### Example: Gram - Schmidt orthogonalization

In this answer, I gave a possible implementation of the Gram - Schmidt orthogonalization procedure, as a higher-order function

GSOrthoNormalizeGen[startvecs_List, dotF_, plusF_, timesF_]

which takes the functions implementing dot product, addition of vectors and multiplication of a vector by a scalar, as parameters. As a result, it can be used for vectors from any vector spaces - all one have to do is to implement these specific functions. In the linked post there are examples for the space of 3D vectors and space of functions.

What this means in practice is that the generic implementation and the specific parts parametrizing it are completely decoupled, they can (and probably should) live in different parts of the project, or even belong to different sub-projects. Therefore, such generalization actually simplifies the code even if I am only interested in a single type of a vector space.

## Closures

Closures are functions that are created at run-time, and have access to the enclosing environment (variables and functions from it). They can then operate on that environment long after the code execution leaves it. Closures are an effective tool to factor and separate functionality. They realize a form on encapsulation, somewhat similar to objects, but more lightweight - they encapsulate behavior rather than state (although can manipulate the state too).

### Example: approximate derivative of a function

This is a classic example. Here, we construct a function that would approximately compute a derivative of another function, numerically:

approxD[f_, dx_] := Function[x, (f[x + dx] - f[x])/dx]

We can now define it for some function:

dsin = approxD[Sin, 0.2];

and use it:

dsin[0]

(* 0.993347 *)

The whole point is that we don't have to know how dsin was constructed, and all the information associated with the process of its construction - we can just use it. It can be stored somewhere, and then used at some later point, perhaps by a different part of the system. It is this kind of separation of construction and execution of certain behavior that makes closures so effective at factoring apart separate components of the system.

### Example: iterator for Fibonacci numbers

Here is a very simple example of a closure, implementing an iterator of Fibonacci numbers:

makeFibIterator[] :=
Module[{current = 1, prev = 0},
Function[Last[{prev, current} = {current, current + prev}]]
]

It is closed over local variables current and prev. It retains access to these variables even after the code leaves Module. This is a stateful closure. Here is an example of use:

iter =  makeFibIterator[];
Table[iter[], 10]

(* {1, 2, 3, 5, 8, 13, 21, 34, 55, 89} *)

Once again, note that after we construct an iterator, it can be used at some point later on, perhaps by a completely different part of the system. That part may not care what this iterator is made of, or how it was constructed - it only knows that the iterator returns a next element when called.

### Summary

Closures are a very useful abstraction, that allows one to encapsulate behaviors and use them later on, without bothering about the full execution context for those behaviors (which may no longer be available), since closures still do have access to that context. They can be thought of as very light-weight objects, but the key difference is that they encapsulate behavior rather than state (but can also manipulate internal state if that's the part of the behavior). They usually work together with higher - order functions, being passed to them and called by them.

Closures facilitate information-hiding, because they allow different parts of the system to exchange minimal units of encapsulated behavior, and the less different parts must assume about the other parts, the more they get decoupled, and the more modular the entire code base becomes. At the same time, they protect their internal state much better than the full-fledged objects in the OO paradigm (unless those have no setters, i.e. are read-only), since they don't provide means to change their internal state other than what they do themselves when being executed. The real need of objects arises when one needs more than one closure with a shared mutable state, but in many cases this isn't really necessary and single closures can do just fine.

Defining abstract data types and making the code de facto more strongly typed is an important tool to scale to larger code bases. It allows to automatically exclude large classes of bugs, which otherwise are likely to appear in the course of code evolution. It also often makes the code much more readable, and much easier to reason about.

There are several possibilities for enforcing stronger typing in one's code. Basically, one can:

• Use patterns and argument checks to enforce types

This is easier to do and is a less formal way to introduce types, which is used in many practical situations. The typing it introduces is similar to a "duck typing", however.

• Use dedicated inert heads as data containers / types

This method allows one to create truly strong types. It might be an overkill do always do that, but there are cases where this option is great.

### Example: using patterns

I will use a function to pick numbers in a specified interval, from this answer

window[list_,{xmin_,xmax_}]:=
Pick[list,Boole[xmin<=#<=xmax]&/@list,1]

We may make it better check arguments, de-facto making this piece of code more strongly-typed:

windowTyped[list_List,{xmin_?NumericQ,xmax_?NumericQ}]:=
Pick[list,Boole[xmin<=#<=xmax]&/@list,1]

We may restrict the arguments even more:

windowTypedStronger[list:{___?NumericQ},{xmin_?NumericQ,xmax_?NumericQ}]:=
Pick[list,Boole[xmin<=#<=xmax]&/@list,1]

Note that you can go further and introduce patterns for your types, so that you can use them in many functions to check the arguments. For example, we could've done:

numericListPtn = {___?NumericQ}

and then

windowTypedStronger[list:numericListPtn, {xmin_?NumericQ,xmax_?NumericQ}]:= ...

where you can now use numericListPtn for argument checks in all functions which are supposed to take this type as an argument

### Example: implementing a Cache data type

This section will illustrate how to implement truly strong types. To make this section a bit more useful, I will post here an implementation of a Cache data structure, based on associations.

Here goes the main implementation:

ClearAll[Cache, RemoveCache, ClearCache];
SetAttributes[Cache, HoldAll];
Cache[limit_]:= Module[{s=<||>},Cache[s, limit]];
Cache /: KeyExistsQ[Cache[s_, _], key_]:=KeyExistsQ[s,key];
Cache /: Cache[s_,_][key_]:=
With[{res = s[key]},
If[!MissingQ[res],AppendTo[s,key->res]];
res
];
Cache /: Normal[Cache[s_,_]]:=s;
(* Note: Append mutates the state and has different semantics here than usual *)
Cache /: Append[c:Cache[s_,limit_], key_-> value_]:=
Module[{},
If[Length[s]==limit,s = Rest[s];];
AppendTo[s,key -> value];
value
];
RemoveCache[Cache[s_,_]]:=CompoundExpression[Remove[s]];
ClearCache[Cache[s_,_]]:= s = <||>;

There are a few things to note from this code:

• The cache objects are strongly typed, they all have the head Cache.
• The Cache head is an inert container for data
• A number of built-in functions were overloaded on Cache data type, using UpValues. In this way, we can use familiar function names without a danger to affect other functionality in the system

The technical part of this implementation is fairly simple. We use the fact that associations are ordered, and when we add a new key-value pair, it is added at the end. The cache is supposed to store n most recent values. To do that, it works as follows: when a value is requested from the cache, and is present there, it moves it at the end (adding the same value again - it is O(1) operation). When we grow the cache to its full capacity, it starts removing key - value pairs from the start. The only tricky part was to have such removal as a fast operation. As nicely pointed out by Mr.Wizard in comments, Rest is O(1), so we use it. Previously, I missed this observation on Rest and used a user-defined analog of Rest here. Note that Delete and Drop on an association are O(n) even for the first positions).

Here are some examples of use:

Create a cache

cache = Cache[10];

Do[Append[cache, i -> i + 1], {i, 1, 10}]

Check what's inside

Normal@cache

(* <|1 -> 2, 2 -> 3, 3 -> 4, 4 -> 5, 5 -> 6, 6 -> 7, 7 -> 8, 8 -> 9, 9 -> 10, 10 -> 11|> *)

Append a key that's already there:

Append[cache, 1 -> 100]

(* 100 *)

We can see that it moved to the right

Normal@cache

(* <|2 -> 3, 3 -> 4, 4 -> 5, 5 -> 6, 6 -> 7, 7 -> 8, 8 -> 9, 9 -> 10, 10 -> 11, 1 -> 100|> *)

Append another key-value

Append[cache, 15 -> 30];
Normal@cache

(* <|3 -> 4, 4 -> 5, 5 -> 6, 6 -> 7, 7 -> 8, 8 -> 9, 9 -> 10, 10 -> 11, 1 -> 100, 15 -> 30|> *)

Do a massive appending (the capacity is only 10)

Do[Append[cache, i -> i + 1], {i, 11, 1000}] // AbsoluteTiming

(* {0.023462, Null} *)

Check the current cache state

Normal@cache

(* <|991 -> 992, 992 -> 993, 993 -> 994, 994 -> 995, 995 -> 996, 996 -> 997, 997 -> 998, 998 -> 999, 999 -> 1000, 1000 -> 1001|> *)

Remove the cache

RemoveCache@cache

### The role of UpValues

It is important to stress that UpValues are an indispensable tool to overload various functions (built-in or user-defined) on custom data types. They provide the only safe way to do that, in fact.

### Summary

Introducing some level of typing is a very useful technique to improve the robustness of your code.

A simpler way to do that is to introduce patterns which expressions belonging to some type should match, and then insert argument checks based on these patterns, into the definitions of those functions which work with these objects. This has an advantage of being simple and quick to do, but a disadvantage that the types introduced in this way won't be in general fully strong and robust.

A somewhat more formal way is to introduce a special head for the type, and then methods working on that head. This is somewhat harder to implement than the first option, but has several advantages: the code is typically more robust and also usually ends up easier to read and understand.

• You will need an ISBN soon. :) And honestly I would buy your book... – gwr Mar 15 '16 at 6:54
• @gwr Thanks :) In fact, I've had plans to write one, may be more than one, for quite some time. Hopefully, those plans will materialize some day. – Leonid Shifrin Mar 15 '16 at 12:38
• I cannot find a case where deleteFirst is faster than Rest; would you give me an example? – Mr.Wizard Mar 16 '16 at 11:29
• @Mr.Wizard Thanks for this observation! For reasons beyond me, I did not try Rest, but only tried Drop and Delete, both of which were slow. I will edit this in. – Leonid Shifrin Mar 16 '16 at 15:55
• Looking again it seems that Delete is also fast in 10.1.0 under Windows, though Drop suffers as described. big = AssociationThread[Range@1*^6, RandomReal[1, 1*^6]]; then Nest[Delete[1], big, 100000] // Length // AbsoluteTiming gives {0.372684, 900000}. The operator form is slightly faster than Delete[#, 1] & but they appear to have similar complexity. – Mr.Wizard Mar 16 '16 at 20:42

This answer is based on the original poster's statement that he has been using Mathematica for six months, and is now trying to build something a bit more complex. I do not take this to mean a large project in the sense that an application developer would use the word.

The notebook interface is really easy to experiment in, and I know that when I had used Mathematica for six months I still needed to test my code very often to make sure it worked. So I used to have a strong preference for the notebook interface. Now I can write much larger pieces of code and be reasonably sure that they will work, so I'm warming up more to IntelliJ which I use sometimes, especially when I am authoring a package.

Unfortunately trying to use both simultaneously - writing on a package in IntelliJ and experimenting with aspects of the package in Mathematica - creates namespace issues that I am usually not willing to fight against, although I have done that too on occasion.

Packages solve mainly three problems in my opinion

1. Code reuse
2. Code distribution
3. Namespace management

If none of these is a concern of yours then you are probably better off writing a notebook in my experience. Note that your mileage with a notebook depends on how you structure it, for larger projects (but not large enough for an IDE) I make use of titles, subtitles and text a lot.

If you find that you have a large code block in a single function then you are probably not writing good, functional code. When you write functional code you want to build up your application gradually using only small functions. The smaller a function is the more testable it is, and this is very important when you are working on a larger project. Otherwise you will find yourself in a situation where you are progressing very slowly because it's hard to debug problems that crop up.

Here is an example of a project with fifty or so functions, but because they have been nested into sensible categories they are all easy to find:

• Nice description/advice. I have followed a similar notebook layout pattern in my development of smaller projects. Also, as mentioned in the post, when I use packages and IntelliJ it is mostly for distribution and reuse purposes. – Anton Antonov Mar 13 '16 at 13:07

As Leonid mentions in one of his answers one of the methods of managing complexity is using Domain Specific Languages (DSLs). In this answer I will provide links to documents, packages, blog posts, and discussions of creating and utilizing DSLs in Mathematica.

For a 2.5 minutes introduction see this video between 25:00 and 27:30.

When to apply DSLs

Here are some situations for applying DSLs.

1. When designing conversational engines.

2. When there are too many usage scenarios and tuning options for the developed algorithms.

• For example, we have a bunch of search, recommendation, and interaction algorithms for a dating site. A different, User Experience Department (UED) designs interactive user interfaces for these algorithms. We make a natural language DSL that invokes the different algorithms according to specified outcomes. With the DSL the different designs produced by UED are much easily prototyped, implemented, or fleshed out. The DSL also gives to UED easier to understand view on the functionalities provided by the algorithms.
3. When designing an API for a collection of algorithms.

Designing DSL

1. Decide what kind of sentences the DSL is going to have.

• Are natural language sentences going to be used?

• Are the language words known beforehand or not?

2. Prepare, create, or accumulate a list of representative sentences.

• In some cases using Morphological Analysis can greatly help for coming up with use cases and the corresponding sentences.
3. Create a context free grammar that describes the sentences from the previous step. (Or a large subset of them.)

4. Program parser(s) for the grammar.

• I use most of the time functional parsers.

• The package FunctionalParsers.m provides a Mathematica implementation of this kind of parsing.

• The package can automatically generate parsers from a grammar given in EBNF. (See the coding example below.)

• I have programmed versions of this package in R and Lua.

5. Program an interpreter for the parsed sentences.

• At this stage the parsed sentences are hooked to the algorithms of the problem domain.

• The package FunctionalParsers.m allows this to be done fairly easy.

6. Test the parsing and interpretation.

See the code example below illustrating steps 3-6.

Introduction to using DSLs in Mathematica

1. This blog post "Natural language processing with functional parsers" gives an introduction to the DSL application in Mathematica.

2. This detailed slide-show presentation "Functional parsers for an integration requests language grammar" shows how to use the package FunctionalParsers.m over a small grammar.

The blog post "Simple time series conversational engine" discusses the creation (design and programming) of a simple conversational engine for time series analysis (data loading, finding outliers and trends.)

Here is a movie demonstrating that conversational engine: http://youtu.be/wlZ5ANglVI4.

Other discussions

1. This answer to the question "How to parse a clojure expression?" provides two concise examples of programming and/or generating parsers for a small grammar.

2. A small part, from 17:30 to 21:00, of the WTC 2012 "Spatial Access Methods and Route Finding" presentation shows a DSL for points of interest queries.

3. The answer of the MSE question "CSS Selectors for Symbolic XML" uses FunctionalParsers.m.

Coding example

This example is for the steps 3-6 of the second section.

 Import["https://raw.githubusercontent.com/antononcube/MathematicaForPrediction/master/FunctionalParsers.m"]

Give an EBNF description of a DSL for food craviings:

ebnfCode = "
<lovefood> = <subject> , <loveverb> , <object-spec> <@ \
LoveFood[Flatten[#]]& ;
<loveverb> = ( 'love' | 'crave' | 'demand' ) <@ LoveType ;
<object-spec> = ( <object-list> | <object> | <objects> | <objects-mult> ) \
<@ LoveObjects[Flatten[{#}]]& ;
<subject> = 'i' | 'we' | 'you' <@ Who ;
<object> = 'sushi' | [ 'a' ] , 'chocolate' | 'milk' | [ 'an' ] , 'ice' , \
'cream' | 'a' , 'tangerine' ;
<objects> = 'sushi' | 'chocolates' | 'milks' | 'ice' , 'creams' | \
'ice-creams' | 'tangerines' ;
<objects-mult> = 'Range[2,100]' , <objects> <@ Mult ;
<object-list> = ( <object> | <objects> | <objects-mult> ) , { 'and' \
\[RightTriangle] ( <object> | <objects> | <objects-mult> ) } ; ";

Generate parses from EBNF string:

GenerateParsersFromEBNF[ToTokens@ebnfCode];

Test the parser pLOVEFOOD for the highest level rule ( ) with a list of sentences:

sentences = {"I love milk", "We demand 2 ice creams",
"I crave 2 ice creams and 5 chocolates",
"You crave chocolate and milk"}; ParsingTestTable[pLOVEFOOD,
ToLowerCase@sentences, "Layout" -> "Horizontal"]

Note the EBNF rule wrappers -- those are symbols specified at the ends of some of the rules.

Next we implement interpreters. I am using WolframAlpha to get the calories. I gave up figuring out how to use EntityValue["Food",___], etc. (Since using WolframAlpha is slow it can be overridden inside Block.)

LoveObjectsCalories[parsed_] :=
Block[{res, wares(*, WolframAlpha={}&*)},
res = (StringJoin @@
Flatten[Riffle[parsed, " and "] /.
Mult[{x_, y_}] :> (StringJoin @@
Riffle[Flatten[{ToString[x], y}], " "])]);
wares = WolframAlpha[res <> " calories", "DataRules"];
{{"Result", 1}, "ComputableData"} /. wares
/. {{"Result", 1}, "ComputableData"} ->
Quantity[RandomInteger[{20, 1200}], "LargeCalories"]
];

LoveFoodCalories[parsed_] :=
Block[{who, type},
who = Cases[parsed, Who[id_] :> id, \[Infinity]][[1]];
type = Cases[parsed, LoveType[id_] :> id, \[Infinity]][[1]];
Which[
who == "you",
Row[{"No, I do not. I am a machine."}],
type == "love",
Row[{"you gain ", Sqrt[1*10.] parsed[[-1]], " per day"}],
True,
Row[{"you will gain ", parsed[[-1]]}]
]
];

Here the parsing tests are done by changing the definitions of the wrapping symbols LoveFood and LoveObjects:

Block[{LoveFood = LoveFoodCalories, LoveObjects = LoveObjectsCalories},
ParsingTestTable[pLOVEFOOD, ToLowerCase@sentences, "Layout" -> "Vertical"]
]

• Thanks a lot, this is just great! I've had it on my todo list to read about parser combinators, but I never had the time, since I'd need to first learn Haskell or ML. And, actually, I've been interested in conversational engines for at least a couple of years, particularly in the context of Mathematica - they seem to have almost endless area of business applications. But I didn't have enough time and didn't know where to start. I had no idea that you published so much great material on the subject! I will certainly find the time to study this. Thanks again! – Leonid Shifrin Mar 16 '16 at 16:19
• @LeonidShifrin Thanks Leonid! At this point I have used functional parsers in many projects. E.g. I have programmed versions of the FunctionalParsers.m package in R and Lua. The Lua implementation allowed me to make DSLs for systems embedded in mobile devices. (E.g. Corona.) I have been looking into denotational semantics and lattice theory in order create and manage families of conversational engines in Mathematica. – Anton Antonov Mar 16 '16 at 17:11
• Sounds very interesting! I wonder what would be the best learning path for someone like me (I only know a little about general principles behind this stuff), to get into this topic, in your opinion. Can you recommend some set of sources (may be books, blog posts, may be some open code / projects, whatever), which one could use to at least get some introductory-level understanding of this (I mean, besides your excellent posts - or do you think your posts could be enough to learn the basics about this)? B.t.w., if that is not something you feel comfortable sharing, I will perfectly understand. – Leonid Shifrin Mar 16 '16 at 17:26
• @LeonidShifrin I need to think about a reading list for someone like you. :) I will share it in a day or two. Some of my friends at WRI have had similar requests. – Anton Antonov Mar 16 '16 at 21:17
• Thanks a bunch! Looking forward! – Leonid Shifrin Mar 16 '16 at 21:19

## Use an IDE like Workbench and remember good software development practices

Take a look at workbench which is a branded version of Eclipse, a very common Integrated Development Platfom (IDE). On the WRI site a lot of information is given of which I would like to point you to a White Paper on Building Large Software in Mathematica. It contains very helpful advice on organizing around packages, using unit tests, version control and a lot more:

Building Large Software in Mathematica (White Paper, 2010)

## "Nothing is more practical than a good theory" (Kurt Lewin)

Many people point at Roman Maeder's book The Mathematica Programmer and I would also add his other book, Computer Science with Mathematica, which I have found to give excellent theoretical background about general principles of programming rather than simply listing some techniques (e.g. know why).

Most helpful there -- and not too commonly found in other books to my knowledge -- are theses topics:

• Structure of Programs (e.g. Packages, Contexts)
• Abstract Data Types (separating algorithmic specifications from its implementation)
• Complexity of Algorithms
• Object-Oriented Programming

(This answer is complementary to those of faysou, C. E., and Leonid Shifrin. It corresponds to one of the methods mentioned but not covered in the answer of Leonid Shifrin.)

For large development projects it is a good idea to use the well established, understood, and documented Object-Oriented Design Patterns as explained and exemplified in my presentation "Object Oriented Design Patterns" at the Wolfram Technology Conference 2015. (The presentation recording is also uploaded at YouTube.)

This diagram shows the large context of patterns:

The talk discusses the use of "Design Patterns by GoF" (in dark red) within Mathematica.

Documented implementation

Here is a link to a document describing how to implement OOP Design Patterns in Mathematica:

"Implementation of Object-Oriented Programming Design Patterns in Mathematica"

The document has almost all of the material covered in my presentation linked above. Of course it has some additional material.

Comparison of Functional Programming and Design Patterns implementations

This MSE answer of mine provides a discussion and links to two implementations of GitHub repository commit history visualization. One of the implementations is with functional programming, the other with Design Patterns. The videos and the document linked above use for examples the package GitHubDataObjects.m from that answer.

• What does the darker green circle to the right mean to someone not having mastered the relevant Asian language? ;) – gwr Mar 13 '16 at 15:28
• @gwr Ha! Good question! It is a reference to the book "Go Rin no Sho" ("The book of five scrolls") by Miyamoto Musashi. We can say that the book contains patterns applicable in antagonistic conflicts. It presents the patterns in abstract forms applicable to person-to-person fights, battles between armies, or other antagonistic settings. Another interesting idea in this book is that if you practice the explained strategy with a samurai sword you will become capable applying that strategy when leading an army. – Anton Antonov Mar 13 '16 at 15:54
• @LeonidShifrin My talk is somewhat rushed and non-linear -- I did not have enough time to cover things in better detail. I hope by end of the day to finish and upload a satisfactory first version of a document describing implementations of OO Design Patterns in Mathematica. (I have worked on it on-and-off for some time...) – Anton Antonov Mar 13 '16 at 18:59
• Well, the Japanese imported Chinese characters for their own use. I'm surprised that the on'yomi and the kun'yomi are the same, as @Shutao notes. – J. M. will be back soon Mar 14 '16 at 2:20
• 我可以閱讀和講一點點, @Shutao. :) – J. M. will be back soon Mar 14 '16 at 2:52

There is a Wolfram Language Tutorial on Wolfram website.

Building Large Software Systems in the Wolfram Language

• Divide the System into Components
• Write and Use Unit Tests
• Think of the Architecture, Not the Code
• Use Wolfram Workbench
• Use Wolfram Language Code Packaging
• Take Advantage of the Wolfram Language
• Keep Things Simple
• Think of Other Developers
• Use Source Control

format: plain


text:
# Linearly independent subsets of a free module

Is it true that the cardinality of every maximal linearly independent subset of a finitely generated free module $A^{n}$ is equal to $n$ (not just at most $n$, but in fact $n$)? Here $A$ is a nonzero commutative ring. I know that it's true if $A$ is Noetherian or integral domain. I thought it was not true in general but I came up with something that looks like a proof and I can't figure out where it went wrong.

-
Dear kwan, I have just deleted an "answer" that didn't address your actual question. My sincere apologies for the misunderstanding, which was entirely my fault. –  Georges Elencwajg Jul 25 '10 at 17:35
I was in fact kicking myself for using the word "rank" and confusing everyone. –  ashpool Jul 25 '10 at 17:42
Georges, I did the same before you did :-) And most other contributors have also fallen into the same trap :-( –  Robin Chapman Jul 25 '10 at 18:08
kwan, I have taken the liberty to remove the word "rank" from the title. I hope that this will help readers discern your true intent. –  Robin Chapman Jul 25 '10 at 18:15
It says the result is false (without giving an example) at this link: books.google.com/… –  KConrad Jul 25 '10 at 20:06

I think I have a counter-example. Let $A$ be the ring of functions $f$ from $\mathbb{C}^2 \setminus (0,0) \to \mathbb{C}$ such there is a polynomial $\widetilde{f} \in \mathbb{C}[x,y]$ such that $\widetilde{f}(x,y)=f(x,y)$ for all but finitely many $(x,y)$ in $\mathbb{C}^2$.

Map $A$ into $A^2$ by $f \mapsto (fx, fy)$. We check that this is injective: If $fx=0$ then $f$ is zero off of the $x$-axis. Similarly, if $fy=0$, then $f$ is zero off of the $y$-axis. So $(fx, fy) = (0,0)$ implies that $f$ is zero everywhere on $\mathbb{C}^2 \setminus (0,0)$.

We now claim that there do not exist $(u,v)$ in $A^2$ such that $(f,g) \mapsto (fx+gu, \ fy+gv)$ is injective. Suppose such a $(u,v)$ exists. Let $\widetilde{u}$ and $\widetilde{v}$ be the polynomials in $\mathbb{C}[x,y]$ which coincide with $u$ and $v$ at all but finitely many points. Let $\Delta=\widetilde{u} y - \widetilde{v} x$. Since $\Delta$ is a polynomial which vanishes at $(0,0)$, it is not a non-zero constant. Thus, $\Delta$ vanishes on an entire infinite subset of $\mathbb{C}^2$. Let $(p,q)$ be a point in $\mathbb{C}^2 \setminus (0,0)$ such that $\Delta(p,q)=0$, $\widetilde{u}(p,q)= u(p,q)$ and $\widetilde{v}(p,q)=v(p,q)$.

So $q u(p,q) - p v(p,q) =0$. Since $(p,q) \neq (0,0)$, there is some $k \in \mathbb{C}$ such that $(u(p,q), v(p,q)) = (kp, kq)$. Take $f$ to be $-k$ at $(p,q)$ and $0$ elsewhere; let $g$ be $1$ at $(p,q)$ and $0$ elsewhere. So $(fx+gu, fy+gv)=0$, and the map $(f,g) \mapsto (fx+gu, \ fy+gv)$ is not injective.

-
Looks good. How did you come up with this? –  KConrad Jul 25 '10 at 20:36
After a lot of dead ends, I decided to concentrate on $(m,n)=(1,2)$. So I wanted a ring with two elements $(x,y)$ such that nothing nonzero is annihilated by the ideal $(x,y)$, but such that every element of the ideal $(x,y)$ annihilates something nonzero. I decided to build $A$ as a subring of $\mathbb{C}^Z$ for some set $Z$. I needed there to be no nonzero function supported on $\{ x=y=0 \}$, but a nonzero function supported on $\{ ay-bx =0 \}$ for any $(a,b) \in A^2$. I tried to find clever ways to do this for a while, before realizing I could just impose it by fiat. –  David Speyer Jul 25 '10 at 20:48
This extends to other examples. Let A be 2-dim. Noeth. domain with inf. many max. ideals. Pick a max. ideal M in A generated by two elements, and let R be the sequences (a_m) indexed by max. ideals m other than M, where a_m is in A/m and there is an a in A such that a_m = a mod m for all but finitely many m (i.e., R is seq. looking like reduction mod m of an elt. of A for all but fin. many m). The ring A embeds in R, so A^2 embeds in R^2. Let M have generators x and y. In R^2, the vector (x,y) is lin. indep. over R but is not part of a 2-elt. lin. indep. subset of R^2. Use C[x,y], Z[x],... –  KConrad Jul 25 '10 at 23:49
Very nice examples! –  Martin Brandenburg Jul 26 '10 at 7:47
Bravo, David. This is a perfect illustration of our site working at its best: ask a rather difficult question and get a clever answer within a few hours. Long live MathOverflow! –  Georges Elencwajg Jul 26 '10 at 19:29

We have to prove that $m \leq n$ if there is a monomorphism $A^m \to A^n$. Since this is given by a $n \times m$ matrix with entries in $A$ and every finitely generated ring is noetherian, it is enough to consider the case that $A$ is noetherian.

Now you already know the proof for this case, but I just add it. Pick a minimal prime ideal $\mathfrak{p} \subseteq A$. This exists since $A \neq 0$. Now localize at $\mathfrak{p}$. Then we may replace $A$ by $A_{\mathfrak{p}}$, and thereby assume that $A$ is a $0$-dimensional noetherian ring, thus artinian. For such a ring it is known that the length of finitely generated modules is finite, and additive on short exact sequences. In particular $m * l(A) \leq n * l(A)$. Since $l(A) \neq 0$ is finite, we get $m \leq n$.

By the way, the assertion can be generalized to the infinite case:

Let $M$ be a free module with basis $B$ and $L \subseteq M$ a linearly independent subset. Then $|L| \leq |B|$.

Proof: Let $B$ be infinite. Representing elements of $L$ as linear combinations of elements in $B$ yields a map $f : L \to E(B)$, where $E(B)$ denotes the set of finite subsets of $B$. Now let $F$ be such a finite subset with $n$ elements. The finite case yields that there are at most $n$ linearly independent elements in $\langle F \rangle$, thus also in $f^{-1}(F)$. Now we use cardinal arithmetics:

$|L| = \sum_{n > 0} \sum_{F \in E(B), |F|=n} |f^{-1}(F)| \leq \sum_{n > 0} |B^n| = \sum_{n > 0} |B| = |B|.$

-
Martin, how does the inequality imply every maximal linearly independent subset have the same cardinality $n$? –  ashpool Jul 25 '10 at 17:54
Ok I realize this is more complicated than I thought ... cf. your recent question mathoverflow.net/questions/30066/…. With the proof above we know that every maximal linearly independent subset of $A^n$ has at most $n$ elements. If it had $m<n$ elements, we should enlarge it somehow ... –  Martin Brandenburg Jul 25 '10 at 18:02
One has to prove more than that an injection $A^m\to A^n$ entails that $m\le n$. One has to show that any linearly independent set of $m<n$ elements in $A^n$ is a subset of a linearly independent set of $n$ elements. –  Robin Chapman Jul 25 '10 at 18:06
Yes, see also my previous comment. Perhaps we can use the Noetherian case ... –  Martin Brandenburg Jul 25 '10 at 18:20
Incidentally the infinite case you mentioned can be also inferred backward from the fact that every maximal linearly independent subset of a free module of infinite rank has the same cardinality. –  ashpool Jul 25 '10 at 22:19

Assuming that $A$ has a maximal ideal $\mathfrak{m}$ (for example, by using Zorn's Lemma), one can proceed as follows: if $M$ is a free $A$-module with basis $(v_i)_{i\in I}$, then $M \cong A^I$, whence $M / \mathfrak{m} M \cong A^I / \mathfrak{m} A^I \cong (A / \mathfrak{m} A)^I$. This is a vector space over $k := A / \mathfrak{m} A$ of dimension $|I|$. Since over fields, all vector space bases of the same vector space have the same length, and since the $k$-vector space structure of $M / \mathfrak{m} M$ is independent of the choice of the basis, this shows that all $A$-bases of $M$ have the same cardinality.

I don't remember where I first saw this though... maybe someone else has a reference? I saw this first in the case that $A = \mathbb{Z}$ and $\mathfrak{m} = (2)$ for free abelian groups $M$, to show that the rank is well-defined.

-
This doesn't show that any linearly independent set of $m<n$ elements in $A^n$ is a subset of a linearly independent set of $n$ elements. –  Robin Chapman Jul 25 '10 at 18:04
Oh, shoot. I misread the question. Thanks for pointing that out, and sorry for my spam... –  felix Jul 25 '10 at 18:12

I have spotted the mistake in my proof. So here is the "wrong" proof:

Let $v_{1},\ldots,v_{m}$ be linearly independent elements in $A^{n}$, where $m\lt n$. Write them as $n$-tuples of elements in $A$, thereby forming an $n$-by-$m$ matrix. Linear independence of $v_{i}\$s means that the rank of this matrix is $m$. So there is an $m$-by-$m$ minor with non-zero determinant. By exchanging rows if necessary, bring these $m$ rows to the top part of the matrix. Now add a colomn to the right side of the matrix whose entries are $0$ except at the $n+1$ th position, where the entry is $1$. Then the new $n$-by-$(m+1)$ matrix has rank $(m+1)$ and hence the $(m+1)$ columns, the first $m$ of which are the $v_{i}\$s, are linearly independent.

The mistake was the notion of rank of a matrix. When the entries are not from integral domain, the proper definition should be the largest integer $m$ such that there is no nonzero element in the ring annihilating the determinant of every $m$-by-$m$ minor. In the above example, I can't conclude that the rank of the $n$-by-$(m+1)$ matrix is $(m+1)$.

With this, I can now exhale a sigh of relief and continue believing that this is not true in general. (By the way I also know that it is true for free modules of infinite rank)

-
continue believe it is not true in general? don't you mean continue believing that it is true in general? Anyway, see my answer, which settles the problem. –  KConrad Jul 25 '10 at 19:14
I have a strong reason to believe that it is not true in general. See mathoverflow.net/questions/30066/… –  ashpool Jul 25 '10 at 19:21
Ah, do you think it is not true because Lazarus did not include "finite free modules" in his list of cases where it is true? –  KConrad Jul 25 '10 at 19:34
Yes. Also see books.google.com/… for non-commutative case. –  ashpool Jul 25 '10 at 19:41