text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Marcos L. Alencastro,  Estefane S. Horn (Revised Gerson B. Bicca)
#* $Id: expander.mso 687 2008-11-20 19:42:33Z bicca $
#*--------------------------------------------------------------------
type expander
	expander()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of chemical components",
				:Lower=>1
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Constant of Gases",
				:Unit=>"kJ/kmol/K",
				:Default=>8.31451,
				:Hidden=>true
			)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Molar Weight",
				:Hidden=>true
			)),(NComp)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Mechanical efficiency",
				:Lower=>1E-3
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Isentropic efficiency",
				:Lower=>1E-3
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Pressure Ratio",
				:Lower=>1E-6,
				:Symbol=>"P_{ratio}"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Drop",
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Decrease",
				:DisplayUnit=>"kPa",
				:Symbol=>"P_{incr}"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Power Losses",
				:Lower=>0
			)),
			energy_mass(Dict{Symbol,Any}(
				:Brief=>"Actual Head",
				:Protected=>true
			)),
			energy_mass(Dict{Symbol,Any}(
				:Brief=>"Isentropic Head",
				:Protected=>true
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Fluid Power",
				:Protected=>true
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Brake Power",
				:Protected=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Isentropic Coefficient",
				:Lower=>0.2,
				:Protected=>true
			)),
			temperature(Dict{Symbol,Any}(
				:Brief=>"Isentropic Temperature",
				:Protected=>true
			)),
			enth_mol(Dict{Symbol,Any}(
				:Brief=>"Enthalpy at constant entropy",
				:Hidden=>true
			)),
			molweight(Dict{Symbol,Any}(
				:Brief=>"Mixture Molar Weight",
				:Hidden=>true
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Mass Density at inlet conditions",
				:Lower=>1E-6,
				:Protected=>true
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Mass Density at outlet conditions",
				:Lower=>1E-6,
				:Protected=>true
			)),
			fraction(Dict{Symbol,Any}(
				:Brief=>"Compressibility factor at inlet",
				:Lower=>1E-3,
				:Protected=>true
			)),
			fraction(Dict{Symbol,Any}(
				:Brief=>"Compressibility factor at outlet",
				:Lower=>1E-3,
				:Protected=>true
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet stream",
				:PosX=>0.14,
				:PosY=>0.0,
				:Symbol=>"_{in}"
			)),
			streamPH(Dict{Symbol,Any}(
				:Brief=>"Outlet stream",
				:PosX=>0.83,
				:PosY=>1,
				:Symbol=>"_{out}"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Work Outlet",
				:PosX=>1,
				:PosY=>0.45,
				:Protected=>true
			)),
			[
				:(Outlet.F = Inlet.F),
				:(Outlet.z = Inlet.z),
				:(Mwm = sum(Mw*Inlet.z)),
				:(Outlet.P = Inlet.P * Pratio),
				:(Outlet.P = Inlet.P - Pdrop),
				:(Outlet.P = Inlet.P - Pdecrease),
				:(rho_in = PP.VapourDensity(Inlet.T, Inlet.P, Inlet.z)),
				:(rho_out= PP.VapourDensity(Outlet.T, Outlet.P, Outlet.z)),
				:(hise = PP.VapourEnthalpy(Tisentropic, Outlet.P, Outlet.z)),
				:(Zfac_in = PP.VapourCompressibilityFactor(Inlet.T,Inlet.P,Inlet.z)),
				:(Zfac_out = PP.VapourCompressibilityFactor(Outlet.T,Outlet.P,Outlet.z)),
				:(Head*Mwm = (Outlet.h-Inlet.h)),
				:(Outlet.T = Tisentropic),
				:((Outlet.h-Inlet.h)= (hise-Inlet.h)*IsentropicEff),
				:(PP.VapourEntropy(Tisentropic, Outlet.P, Outlet.z) = PP.VapourEntropy(Inlet.T, Inlet.P, Inlet.z)),
				:(BrakePower = WorkOut),
				:(BrakePower = FluidPower*MechanicalEff),
				:(PowerLoss = BrakePower - FluidPower),
				:(FluidPower = HeadIsentropic*Mwm*Inlet.F*IsentropicEff),
				:(HeadIsentropic*Mwm*((IseCoeff-1.001)/IseCoeff) = (0.5*Zfac_in+0.5*Zfac_out)*Rgas*Inlet.T*((Pratio)^((IseCoeff-1.001)/IseCoeff) - 1)),
				:(HeadIsentropic*Mwm = (hise -Inlet.h)),
			],
			[
				"Overall Molar Balance","Component Molar Balance","Average Molecular Weight","Pressure Ratio","Pressure Drop","Pressure Decrease","Mass Density at inlet conditions","Mass Density at outlet conditions","Enthalpy at isentropic conditions","Compressibility factor at Inlet Conditions","Compressibility factor at Outlet Conditions","Actual Head","Discharge Temperature","Discharge Temperature","Isentropic Outlet Temperature","Brake Power","Brake Power","Power Loss","Fluid Power","Isentropic Coefficient","Isentropic Head",
			],
			[:PP,:NComp,:Rgas,:Mw,:MechanicalEff,:IsentropicEff,],
			[:Pratio,:Pdrop,:Pdecrease,:PowerLoss,:Head,:HeadIsentropic,:FluidPower,:BrakePower,:IseCoeff,:Tisentropic,:hise,:Mwm,:rho_in,:rho_out,:Zfac_in,:Zfac_out,:Inlet,:Outlet,:WorkOut,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Rgas::positive
	Mw::Array{molweight}
	MechanicalEff::positive
	IsentropicEff::positive
	Pratio::positive
	Pdrop::press_delta
	Pdecrease::press_delta
	PowerLoss::power
	Head::energy_mass
	HeadIsentropic::energy_mass
	FluidPower::power
	BrakePower::power
	IseCoeff::positive
	Tisentropic::temperature
	hise::enth_mol
	Mwm::molweight
	rho_in::dens_mass
	rho_out::dens_mass
	Zfac_in::fraction
	Zfac_out::fraction
	Inlet::stream
	Outlet::streamPH
	WorkOut::power
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export expander
function set(in::expander)
	Mw = PP.MolecularWeight()
	 Rgas = 8.31451*"kJ/kmol/K"

end
function setEquationFlow(in::expander)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	if IsentropicEff >= 1
		addEquation(13)
	else
		addEquation(14)
	end
	addEquation(15)
	addEquation(16)
	addEquation(17)
	addEquation(18)
	addEquation(19)
	addEquation(20)
	addEquation(21)
end
function atributes(in::expander,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/expander"
	fields[:Brief]="Model of an expander."
	fields[:Info]="To be documented

== References ==

[1] GPSA, 1979, Engineering Data Book, Chapter 4, 5-9 - 5-10.

[2] Bloch, Heinz P., A Practical Guide to Compressor Technology, John Wiley & Sons, Incorporate, 2006."
	drive!(fields,_)
	return fields
end
expander(_::Dict{Symbol,Any})=begin
	newModel=expander()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 77311686-7a2b-11eb-2bc7-95b2a2211969
begin
	using Pkg
	Pkg.activate("MLJ_env", shared=true)

	# using MLJ
	using CSV
	using DataFrames
	using PlutoUI
	using Random
	using Test
	using Printf
	# using Plots
end

# ╔═╡ 3578d030-7eef-11eb-1d64-d7e7ff02159c
begin
	const DF = AbstractDataFrame # Union{DataFrame, SubDataFrame}

	include("./utils.jl");
	include("./dt_utils.jl");
end

# ╔═╡ 5435ff66-7a2b-11eb-1474-b96dcad21315
md"""
## C03w03: Implementing binary Decision Trees

The goal of this notebook is to implement your own binary decision tree classifier. We will:

  - Use DataFrames to do some feature engineering.
  - Transform categorical variables into binary variables.
  - Write a function to compute the number of misclassified examples in an intermediate node.
  - Write a function to find the best feature to split on.
  - Build a binary decision tree from scratch.
  - Make predictions using the decision tree.
  - Evaluate the accuracy of the decision tree.
  - Visualize the decision at the root node.

**Important Note**: In this assignment, we will focus on building decision trees where the data contain **only binary (0 or 1) features**. This allows us to avoid dealing with:
  -  Multiple intermediate nodes in a split
  - The thresholding issues of real-valued features.
"""

# ╔═╡ 771130dc-7a2b-11eb-170f-934c11fceede
md"""
#### Load LendingClub dataset
"""

# ╔═╡ 76f62e40-7a2b-11eb-095e-f36e207f06a2
begin
	loans =  CSV.File(("../../ML_UW_Spec/C03/data/lending-club-data.csv");
	header=true) |> DataFrame;
	size(loans)
end

# ╔═╡ 76de16d4-7a2b-11eb-09c6-e36b6c893c1f
md"""
Like the previous assignment, we reassign the labels to have +1 for a safe loan, and -1 for a risky (bad) loan.
"""

# ╔═╡ 76c17c4a-7a2b-11eb-2c47-97106bd58f99
insertcols!(loans, :safe_loans => ifelse.(loans.bad_loans .== 0, 1, -1),
 	makeunique=true);

# ╔═╡ 76a83684-7a2b-11eb-1961-cb133db8c164
md"""
Unlike the previous assignment where we used several features, in this assignment, we will just be using 4 categorical
features:

  1. grade of the loan
  2. the length of the loan term
  3. the home ownership status: own, mortgage, rent
  4. number of years of employment.

Since we are building a binary decision tree, we will have to convert these categorical features to a binary representation in a subsequent section using 1-hot encoding.
"""

# ╔═╡ 768d6b08-7a2b-11eb-388a-8763c8da0a51
begin
  	const Features = [
		:grade,          # grade of the loan
		:emp_length,     # number of years of employment
		:home_ownership, # home_ownership status: own, mortgage or rent
		:term,           # the term of the loan
	]

	const Target = :safe_loans # prediction target (y) (+1 means safe, -1 is risky)

	# Extract the feature columns and target column
	select!(loans, [Features..., Target]);
	length(names(loans)), names(loans)
end

# ╔═╡ dde7e8e4-7a2d-11eb-2527-3df1ab7d0796
describe(loans, :eltype, :nmissing, :first => first)

# ╔═╡ 7673d8e4-7a2b-11eb-02b4-2f837181b4e9
md"""
#### Subsample dataset to make sure classes are balanced

Just as we did in the previous assignment, we will undersample the larger class (safe loans) in order to balance out our dataset. This means we are throwing away many data points.
"""

# ╔═╡ 542db682-7a2e-11eb-35ad-25c82d6919e0
begin
	len_loans₀	 = size(loans)[1]
  	safe_loans₀, risky_loans₀ = df_partition(loans, Target)
  	p_safe_loans₀, p_risky_loans₀ = df_ratios(safe_loans₀, risky_loans₀, len_loans₀);

  	with_terminal() do
    	@printf("Number of safe loans  : %d\n", size(safe_loans₀)[1])
    	@printf("Number of risky loans : %d\n", size(risky_loans₀)[1])
    	@printf("Percentage of safe loans:  %1.2f%%\n", p_safe_loans₀ * 100.)
    	@printf("Percentage of risky loans: %3.2f%%\n", p_risky_loans₀ * 100.)
  	end
end

# ╔═╡ c5038900-7eec-11eb-34ff-9b3613ca036c
loans_data = resample!(safe_loans₀, risky_loans₀);

# ╔═╡ 3d3614d8-7a2e-11eb-2a8a-b53822e7928e
begin
	len_loans₁	 = size(loans_data)[1]
  	safe_loans₁, risky_loans₁ = df_partition(loans_data, Target)
  	p_safe_loans₁, p_risky_loans₁ = df_ratios(safe_loans₁, risky_loans₁, len_loans₁);

 	with_terminal() do
    	@printf("Number of safe loans  : %d\n", size(safe_loans₁)[1])
    	@printf("Number of risky loans : %d\n", size(risky_loans₁)[1])
    	@printf("Percentage of safe loans:  %1.2f%%\n", p_safe_loans₁ * 100.)
    	@printf("Percentage of risky loans: %3.2f%%\n", p_risky_loans₁ * 100.)
  	end
end

# ╔═╡ 760c171a-7a2b-11eb-1341-19ba5a1824a1
md"""
#### Transform categorical data into binary features

In this assignment, we will implement **binary decision trees** (decision trees for binary features, a specific case of categorical variables taking on two values, e.g., true/false). Since all of our features are currently categorical features, we want to turn them into binary features.

For instance, the `home_ownership` feature represents the home ownership status of the loanee, which is either `own`, `mortgage` or `rent`. For example, if a data point has the feature
```
   {'home_ownership': 'RENT'}
```
we want to turn this into three features:
```
 {
   'home_ownership = OWN'      : 0,
   'home_ownership = MORTGAGE' : 0,
   'home_ownership = RENT'     : 1
 }
```
"""

# ╔═╡ 75d98df2-7a2b-11eb-296e-ef4d312ffcfa
loans_data[1:5, :]

# ╔═╡ 75c1a64e-7a2b-11eb-21f6-55dd3a8ce164
begin
	hot_encode!(loans_data);
	loans_data[1:5, :];
end

# ╔═╡ 75a52c44-7a2b-11eb-2f02-2312f3a3dac9
begin
	features = names(loans_data) |> a_ -> Symbol.(a_);
	deleteat!(features,
		findall(x -> x == Target, features))
end

# ╔═╡ 44158e34-7a40-11eb-1d64-d7e7ff02159c
first(loans_data, 3)

# ╔═╡ 75736a9c-7a2b-11eb-2278-ab37447cf22d
## Expect 25 features after hot-encoding
@test length(features) == 25

# ╔═╡ 116eb5f8-7a33-11eb-2c4c-5f937540f987
length(loans_data.grade_A)
# vs 46508 in original notebook

# ╔═╡ a5c56d68-7a32-11eb-2ef6-e58ed65dc561
md"""
**Checkpoint:** Make sure the following answers match up.
"""

# ╔═╡ a5a1cdf6-7a32-11eb-18dd-0dfee8eb885c
@test sum(loans_data.grade_A) == 6537
# vs 6422 in original notebook

# ╔═╡ a56e44e0-7a32-11eb-2587-6f5f8b8c21f3
md"""
#### Train-test split

We split the data into a train test split with 80% of the data in the training set and 20% of the data in the test set.
"""

# ╔═╡ 6fad8afe-7a33-11eb-3208-6b289c3ccd23
begin
	train_data, test_data = train_test_split(loans_data; split=0.8, seed=42);
	size(train_data), size(test_data)
    ## vs ((37224, 26), (9284, 26)) ...
end

# ╔═╡ 6f88e4f6-7a33-11eb-10cf-e3e8608e616b
md"""
### Decision tree implementation

In this section, we will implement binary decision trees from scratch. There are several steps involved in building a decision tree. For that reason, we have split the entire assignment into several sections.

#### Function to count number of mistakes while predicting majority class

Recall from the lecture that prediction at an intermediate node works by predicting the **majority class** for all data points that belong to this node.

Now, we will write a function that calculates the number of **missclassified examples** when predicting the **majority class**. This will be used to help determine which feature is the best to split on at a given node of the tree.

**Note**: Keep in mind that in order to compute the number of mistakes for a majority classifier, we only need the label (y values) of the data points in the node.

**Steps to follow**:
  1. Calculate the number of safe loans and risky loans.
  1. Since we are assuming majority class prediction, all the data points that are **not** in the majority class are considered **mistakes**.
  1. Return the number of **mistakes**.

Now, let us write the function `inter_node_num_mistakes` which computes the number of misclassified examples of an intermediate node given the set of labels (y values) of the data points contained in the node.
"""

# ╔═╡ 6f6802b8-7a33-11eb-0c47-5d8b41dcfb39
function inter_node_num_mistakes(labels_in_node::AbstractVector{T}) where {T <: Real}
    length(labels_in_node) == 0 && return 0  ## Corner case

	n = length(labels_in_node)
	n_pos = sum(labels_in_node .== 1)  ## Count the number of 1's (safe loans)
    n_neg = n - n_pos

	## Return the number of mistakes that the majority classifier makes.
	return n_pos ≥ n_neg ? n_neg : n_pos
end

# ╔═╡ 6f4a5e40-7a33-11eb-038f-438ba7d14f44
md"""
Because there are several steps in this assignment, we have introduced some stopping points where you can check your code and make sure it is correct before proceeding.
To test your inter_node_num_mistakes function, run the following code until you get a Test passed!, then you should proceed. Otherwise, you should spend some time figuring out where things went wrong.
"""

# ╔═╡ 6f3011da-7a33-11eb-2cd0-b7bb04b7b545
with_terminal() do
	@testset "begin sanity checks" begin
		## Test case 1
		ex1_labels = DataFrame(:x => [-1, -1, 1, 1, 1])
		@test inter_node_num_mistakes(ex1_labels.x) == 2

		## Test case 2
		ex2_labels = [-1, -1, 1, 1, -1, 1, 1]
		@test inter_node_num_mistakes(ex2_labels) == 3

		## Test case 3
		ex3_labels = DataFrame(:x => [-1, -1, -1, -1, 1, -1, 1])
		@test inter_node_num_mistakes(ex3_labels.x) == 2
	end
end

# ╔═╡ 6f177b72-7a33-11eb-0972-f35a77f51499
md"""
#### Function to pick best feature to split on

The function `best_splitting_feature` takes 3 arguments:
1. The data (DataFrame which includes all of the feature columns and label column)
2. The features to consider for splits (a list of Symbol of column names to consider for splits)
3. The name of the target/label column (Symbol)

The function will loop through the list of possible features, and consider splitting on each of them. It will calculate the classification error of each split and return the feature that had the smallest classification error when split on.

Recall that the **classification error** is defined as follows:

$$\mbox{classification error} = \frac{\mbox{\# mistakes}}{\mbox{\# total examples}}$$

Follow these steps:
  - **Step 1:** Loop over each feature in the feature list
  -  **Step 2:** Within the loop, split the data into two groups: one group where all of the data has feature value 0 or False (we will call this the **left** split), and one group where all of the data has feature value 1 or True (we will call this the **right** split). Make sure the **left** split corresponds with 0 and the **right** split corresponds with 1 to ensure your implementation fits with our implementation of the tree building process.
  - **Step 3:** Calculate the number of misclassified examples in both groups of data and use the above formula to compute the **classification error**.
  - **Step 4:** If the computed error is smaller than the best error found so far, store this **feature and its error**.

This may seem like a lot, but we have provided pseudocode in the comments in order to help you implement the function correctly.

**Note:** Remember that since we are only dealing with binary features, we do not have to consider thresholds for real-valued features. This makes the implementation of this function much easier.
"""

# ╔═╡ 6efa48f4-7a33-11eb-1dac-3bf64871bf99
function best_splitting_feature(df::DF, features::Vector{Symbol},
		target::Symbol) where {DF <: Union{DataFrame, SubDataFrame}}
	best_feature, best_error = nothing, 2
	## Note: Since error is always ≤ 1, init best_error with something > 1

	num_data_points = size(df)[1]
	for f ∈ features
		l_split = df[df[!, f] .== 0, :]
		r_split = df[df[!, f] .== 1, :]

		l_mistakes = inter_node_num_mistakes(l_split[!, target])
		r_mistakes = inter_node_num_mistakes(r_split[!, target])

		error = (l_mistakes + r_mistakes) / num_data_points
		if error < best_error
			best_error = error
          	best_feature = f
		end
	end

	best_feature
end

# ╔═╡ 57fb569a-7a36-11eb-1d64-d7e7ff02159c
@test best_splitting_feature(train_data, features, Target) == Symbol("term_ 60 months")

# ╔═╡ 6ee0a5a2-7a33-11eb-31db-d164581fd26a
md"""
#### Building the tree

With the above functions implemented correctly, we are now ready to build our decision tree. Each node in the decision tree is represented as a dictionary which contains the following keys and possible values:

```julia
    (
       :is_leaf            => true/false.
       :prediction         => Prediction at the leaf node.
       :left               => (dictionary corresponding to the left tree).
       :right              => (dictionary corresponding to the right tree).
       :splitting_feature  => The feature that this node splits on.
    )
```

First, we will write a function that creates a leaf node given a set of target values.
"""

# ╔═╡ a553b5ee-7a32-11eb-329a-419005d73e2c
function create_leaf(target_values;
		splitting_feature=nothing,
		left=nothing, right=nothing, is_leaf=true)

    ## Create a leaf node
    leaf = Dict{Symbol, Any}(
      :splitting_feature => splitting_feature,
	  :prediction => nothing,
      :left => left,
      :right => right,
      :is_leaf => is_leaf,
    )

	if is_leaf
    	## Count the number of data points that are +1 and -1 in this node.
    	num_1s = length(target_values[target_values .== 1, :])
    	num_minus_1s = length(target_values[target_values .== -1, :])

    	## For the leaf node, set the prediction to be the majority class.
    	## Store the predicted class (1 or -1) in leaf['prediction']
    	leaf[:prediction] = num_1s > num_minus_1s ? 1 : -1
	end

	leaf
end

# ╔═╡ b715aca6-7a37-11eb-30dd-677a45cd9920
md"""
We have provided a function that learns the decision tree recursively and implements 3 stopping conditions:
  1. **Stopping condition 1:** All data points in a node are from the same class.
  2. **Stopping condition 2:** No more features to split on.
  3. **Additional stopping condition:** In addition to the above two stopping conditions covered in lecture, in this assignment we will also consider a stopping condition based on the **max_depth** of the tree. By not letting the tree grow too deep, we will save computational effort in the learning process.

Now, we will write down the skeleton of the learning algorithm.
"""

# ╔═╡ b6fc15ae-7a37-11eb-071c-0fb7dd06e885
function decision_tree_create(df::DF, features, target;
		curr_depth=0, max_depth=10,
		verbose=false) where {DF <: Union{DataFrame, SubDataFrame}}

    rem_features = copy(features) # Make a copy of the features.

    target_vals = df[!, target]
    if verbose
      println("-------------------------------------------------------------")
      @printf("Subtree, depth = %s (%s data points).\n", curr_depth,
			size(target_vals)[1])
	end

    if inter_node_num_mistakes(target_vals) == 0
		## Stopping cond1: no mistakes at current node?
        verbose && @printf("Stopping condition 1 reached.")
        return create_leaf(target_vals)

	elseif length(rem_features) == 0
    	## Stopping cond2: no remaining features to consider splitting on?
        verbose && @printf("Stopping condition 2 reached.")
        return create_leaf(target_vals)

	elseif curr_depth ≥ max_depth
    	## Additional stopping condition (limit tree depth)
        verbose && @printf("Reached maximum depth. Stopping for now.")
        return create_leaf(target_vals)
	end

    ## Now Find the best splitting feature
    splitting_feature = best_splitting_feature(df, features, target)

    # Split on the best feature that we found.
    l_split = df[df[!, splitting_feature] .== 0, :]
    r_split = df[df[!, splitting_feature] .== 1, :]

	deleteat!(rem_features,
		findall(x -> x == splitting_feature, rem_features))

    verbose && @printf("Split on feature %s. (%s, %s)\n", splitting_feature,
		size(l_split)[1], size(r_split)[1])

    # Create a leaf node if the split is "perfect"
    if size(l_split)[1] == size(df)[1]
        verbose && println("Creating leaf node.")
        return create_leaf(l_split[!, target])
	end

    if size(r_split)[1] == size(df)[1]
        verbose && println("Creating leaf node.")
        return create_leaf(r_split[!, target])
	end

    # Recurse on left and right subtrees
    ltree = decision_tree_create(l_split, rem_features, target;
						curr_depth=curr_depth + 1,
						max_depth, verbose)
    rtree = decision_tree_create(r_split, rem_features, target;
						curr_depth=curr_depth + 1,
						max_depth, verbose)

    return create_leaf(target_vals;
						splitting_feature=string(splitting_feature),
						left=ltree,
						right=rtree,
						is_leaf=false)
end

# ╔═╡ ae0b1116-7a3a-11eb-37bc-7df8411fb645
# function count_nodes(tree)
#     isnothing(tree) || tree[:is_leaf] ?
# 	  1 : 1 + count_nodes(tree[:left]) + count_nodes(tree[:right])
# end

# ╔═╡ b6e1567c-7a37-11eb-2890-7577dda5d4de
md"""

Run the following test code to check the implementation.
"""

# ╔═╡ b6c80f78-7a37-11eb-34ff-9b3613ca036c
begin
	small_data_decision_tree = decision_tree_create(train_data, features, :safe_loans;
		max_depth=3, verbose=true)

	@test count_nodes(small_data_decision_tree) == 13
end

# ╔═╡ 753b3df2-7a2b-11eb-28c0-b5bdffef414b
md"""
### Build the tree!

Now that all the tests are passing, we will train a tree model on the train_data. Limit the depth to 6 `max_depth = 6` to make sure the algorithm doesn't run for too long.

Call this tree `dt6_model`.

"""

# ╔═╡ be0730f6-7a3d-11eb-30dd-677a45cd9920
dt6_model = decision_tree_create(train_data, features, :safe_loans;
		max_depth=6)

# ╔═╡ bdeea6bc-7a3d-11eb-071c-0fb7dd06e885
md"""
### Making predictions with a decision tree

As discussed in the lecture, we can make predictions from the decision tree with a simple recursive function. Below, we call this function `classify` (cf. `dt_utils.jl`), which takes in a learned tree and a test point x to classify. We include an option annotate that describes the prediction path when set to True.

"""

# ╔═╡ bdb7d658-7a3d-11eb-34ff-9b3613ca036c
md"""
Now, let's consider the first example of the test set and see what `dt6_model`  predicts for this data point.
"""

# ╔═╡ e82bb3ce-7a3e-11eb-1d64-d7e7ff02159c
with_terminal() do
	@printf("Predicted class: %s\n", classify(dt6_model, test_data[1, :];
			annotate=true))
end

# ╔═╡ 96aa449a-7a4b-11eb-37bc-7df8411fb645
md"""
**Quiz Question**: What was the feature that dt6_model first split on while making the prediction for `test_data[1, :]`?
  - "term 60 months"

**Quiz Question**: What was the first feature that lead to a right split of `test_data[1, :]`?
  - "gradeF" (where split\_feature\_value was 1)

**Quiz Question**: What was the last feature split on before reaching a leaf node for `test_data[1, :]`?
  - "gradeF"

"""

# ╔═╡ 94f1a986-7a4b-11eb-30dd-677a45cd9920
md"""
#### Evaluating your decision tree

Now, we will write a function to evaluate a decision tree by computing the classification error of the tree on the given dataset.

Again, recall that the **classification error** is defined as follows:

$$\mbox{classification error} = \frac{\mbox{\# mistakes}}{\mbox{\# total examples}}$$

Now, write a function called `evaluate_classification_error` that takes in as input:
1. `tree` (as described above)
2. `data` (an SFrame)
3. `target` (a string - the name of the target/label column)

This function should calculate a prediction (class label) for each row in `data` using the decision `tree` and return the classification error computed using the above formula.
"""

# ╔═╡ 7b7e1d1c-7a4c-11eb-0a51-851d7325ad47
function eval_classification_error(tree, df::DF, target::Symbol) where {DF <: Union{DataFrame, SubDataFrame}}
	## 1. get the predictions
    ŷ = map(x -> classify(tree, x), eachrow(df))

    ## 2. calculate the classification error
    num_mistakes = sum(ŷ .≠ df[!, target])
    return num_mistakes / size(df)[1]
end

# ╔═╡ 7b66e926-7a4c-11eb-08e3-635d3737c018
begin
	error = eval_classification_error(dt6_model, test_data, Target)
	(class_error=round(error, digits=2),
		class_accuracy=round(1. - error, digits=2))
end

# ╔═╡ 7b4d4fb6-7a4c-11eb-0a69-8fbbd0859786
md"""
**Quiz Question**: Rounded to 2nd decimal point, what is the classification error of `dt6_model` on the test\_data?
  - cf. above cell.
"""

# ╔═╡ 7b31c35e-7a4c-11eb-05df-8b16a1439d68
md"""
#### Printing out a decision stump

As discussed in the lecture, we can print out a single decision stump .
"""

# ╔═╡ 881bdbb2-7a57-11eb-20d9-87257af20d51
function  print_stump(tree, name=:root)
    split_name = tree[:splitting_feature] # ex. 'term. 36 months'

    if isnothing(split_name)
        print("(leaf, label: $(tree[:prediction]))")
        return nothing
	end
	# println("==> ", split_name)
    # split_feature, split_value = split(split_name, "_")
    println("                       $(name)")
    println("         |---------------|----------------|")
    println("         |                                |")
    println("         |                                |")
    println("         |                                |")
    println("     [$(split_name) == 0]              [$(split_name) == 1]")
    println("         |                                |")
    println("         |                                |")
    println("         |                                |")
    @printf("     (%s)                     (%s)\n",
        tree[:left][:is_leaf] ?
			string("leaf, label: ", string(tree[:left][:prediction])) : "subtree",
		tree[:right][:is_leaf] ?
        	string("leaf, label: ", string(tree[:right][:prediction])) : "subtree")
end

# ╔═╡ caea5eb8-7a58-11eb-12fa-4d42c67876b2
with_terminal() do
	print_stump(dt6_model)
end

# ╔═╡ 88b9603c-7a5a-11eb-2ad1-017f6ac940dd
with_terminal() do
	dt6_model[:splitting_feature]
end

# ╔═╡ aca01ac2-7a5a-11eb-072e-f31cc84a982a
with_terminal() do
	print_stump(dt6_model[:left], dt6_model[:splitting_feature])
end

# ╔═╡ ac83442e-7a5a-11eb-12ec-a5cf5ea3c4c1
with_terminal() do
	print_stump(dt6_model[:left][:left][:left],
		dt6_model[:left][:left][:splitting_feature])
end

# ╔═╡ ac6589fc-7a5a-11eb-1bbb-e39f6ba35768
with_terminal() do
	print_stump(dt6_model[:right][:right], dt6_model[:right][:splitting_feature])
end

# ╔═╡ Cell order:
# ╟─5435ff66-7a2b-11eb-1474-b96dcad21315
# ╠═77311686-7a2b-11eb-2bc7-95b2a2211969
# ╠═3578d030-7eef-11eb-1d64-d7e7ff02159c
# ╟─771130dc-7a2b-11eb-170f-934c11fceede
# ╠═76f62e40-7a2b-11eb-095e-f36e207f06a2
# ╟─76de16d4-7a2b-11eb-09c6-e36b6c893c1f
# ╠═76c17c4a-7a2b-11eb-2c47-97106bd58f99
# ╟─76a83684-7a2b-11eb-1961-cb133db8c164
# ╠═768d6b08-7a2b-11eb-388a-8763c8da0a51
# ╠═dde7e8e4-7a2d-11eb-2527-3df1ab7d0796
# ╟─7673d8e4-7a2b-11eb-02b4-2f837181b4e9
# ╠═542db682-7a2e-11eb-35ad-25c82d6919e0
# ╠═c5038900-7eec-11eb-34ff-9b3613ca036c
# ╠═3d3614d8-7a2e-11eb-2a8a-b53822e7928e
# ╟─760c171a-7a2b-11eb-1341-19ba5a1824a1
# ╠═75d98df2-7a2b-11eb-296e-ef4d312ffcfa
# ╠═75c1a64e-7a2b-11eb-21f6-55dd3a8ce164
# ╠═75a52c44-7a2b-11eb-2f02-2312f3a3dac9
# ╠═44158e34-7a40-11eb-1d64-d7e7ff02159c
# ╠═75736a9c-7a2b-11eb-2278-ab37447cf22d
# ╠═116eb5f8-7a33-11eb-2c4c-5f937540f987
# ╟─a5c56d68-7a32-11eb-2ef6-e58ed65dc561
# ╠═a5a1cdf6-7a32-11eb-18dd-0dfee8eb885c
# ╟─a56e44e0-7a32-11eb-2587-6f5f8b8c21f3
# ╠═6fad8afe-7a33-11eb-3208-6b289c3ccd23
# ╟─6f88e4f6-7a33-11eb-10cf-e3e8608e616b
# ╠═6f6802b8-7a33-11eb-0c47-5d8b41dcfb39
# ╟─6f4a5e40-7a33-11eb-038f-438ba7d14f44
# ╠═6f3011da-7a33-11eb-2cd0-b7bb04b7b545
# ╟─6f177b72-7a33-11eb-0972-f35a77f51499
# ╠═6efa48f4-7a33-11eb-1dac-3bf64871bf99
# ╠═57fb569a-7a36-11eb-1d64-d7e7ff02159c
# ╟─6ee0a5a2-7a33-11eb-31db-d164581fd26a
# ╠═a553b5ee-7a32-11eb-329a-419005d73e2c
# ╟─b715aca6-7a37-11eb-30dd-677a45cd9920
# ╠═b6fc15ae-7a37-11eb-071c-0fb7dd06e885
# ╠═ae0b1116-7a3a-11eb-37bc-7df8411fb645
# ╟─b6e1567c-7a37-11eb-2890-7577dda5d4de
# ╠═b6c80f78-7a37-11eb-34ff-9b3613ca036c
# ╟─753b3df2-7a2b-11eb-28c0-b5bdffef414b
# ╠═be0730f6-7a3d-11eb-30dd-677a45cd9920
# ╟─bdeea6bc-7a3d-11eb-071c-0fb7dd06e885
# ╟─bdb7d658-7a3d-11eb-34ff-9b3613ca036c
# ╠═e82bb3ce-7a3e-11eb-1d64-d7e7ff02159c
# ╟─96aa449a-7a4b-11eb-37bc-7df8411fb645
# ╟─94f1a986-7a4b-11eb-30dd-677a45cd9920
# ╠═7b7e1d1c-7a4c-11eb-0a51-851d7325ad47
# ╠═7b66e926-7a4c-11eb-08e3-635d3737c018
# ╟─7b4d4fb6-7a4c-11eb-0a69-8fbbd0859786
# ╟─7b31c35e-7a4c-11eb-05df-8b16a1439d68
# ╠═881bdbb2-7a57-11eb-20d9-87257af20d51
# ╠═caea5eb8-7a58-11eb-12fa-4d42c67876b2
# ╠═88b9603c-7a5a-11eb-2ad1-017f6ac940dd
# ╠═aca01ac2-7a5a-11eb-072e-f31cc84a982a
# ╠═ac83442e-7a5a-11eb-12ec-a5cf5ea3c4c1
# ╠═ac6589fc-7a5a-11eb-1bbb-e39f6ba35768

ext: jl


text:
using Base.Test
using ShowSet

A = IntSet(1,3,9)
@test string(A) == "{1,3,9}"

B = Set{ASCIIString}()
push!(B,"Hello")
push!(B,"Bye")
@test "{Bye,Hello}" == string(B)


text:
using IterativeSolvers
using LinearMaps
using Test
using LinearAlgebra
using SparseArrays
using Random

import LinearAlgebra.ldiv!

include("laplace_matrix.jl")

struct JacobiPrec
    diagonal
end

ldiv!(y, P::JacobiPrec, x) = y .= x ./ P.diagonal

@testset "Conjugate Gradients" begin

Random.seed!(1234321)

@testset "Small full system" begin
    n = 10

    @testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
        A = rand(T, n, n)
        A = A' * A + I
        b = rand(T, n)
        tol = √eps(real(T))

        x,ch = cg(A, b; tol=tol, maxiter=2n, log=true)
        @test isa(ch, ConvergenceHistory)
        @test norm(A*x - b) / norm(b) ≤ tol
        @test ch.isconverged

        # If you start from the exact solution, you should converge immediately
        x,ch = cg!(A \ b, A, b; tol=10tol, log=true)
        @test niters(ch) ≤ 1
        @test nprods(ch) ≤ 2

        # Test with cholfact should converge immediately
        F = cholesky(A, Val(false))
        x,ch = cg(A, b; Pl=F, log=true)
        @test niters(ch) ≤ 2
        @test nprods(ch) ≤ 2

        # All-zeros rhs should give all-zeros lhs
        x0 = cg(A, zeros(T, n))
        @test x0 == zeros(T, n)
    end
end

@testset "Sparse Laplacian" begin
    A = laplace_matrix(Float64, 10, 2)
    P = JacobiPrec(diag(A))

    rhs = randn(size(A, 2))
    rmul!(rhs, inv(norm(rhs)))
    tol = 1e-5

    @testset "Matrix" begin
        xCG = cg(A, rhs; tol=tol, maxiter=100)
        xJAC = cg(A, rhs; Pl=P, tol=tol, maxiter=100)
        @test norm(A * xCG - rhs) ≤ tol
        @test norm(A * xJAC - rhs) ≤ tol
    end

    Af = LinearMap(A)
    @testset "Function" begin
        xCG = cg(Af, rhs; tol=tol, maxiter=100)
        xJAC = cg(Af, rhs; Pl=P, tol=tol, maxiter=100)
        @test norm(A * xCG - rhs) ≤ tol
        @test norm(A * xJAC - rhs) ≤ tol
    end

    @testset "Function with specified starting guess" begin
        x0 = randn(size(rhs))
        xCG, hCG = cg!(copy(x0), Af, rhs; tol=tol, maxiter=100, log=true)
        xJAC, hJAC = cg!(copy(x0), Af, rhs; Pl=P, tol=tol, maxiter=100, log=true)
        @test norm(A * xCG - rhs) ≤ tol
        @test norm(A * xJAC - rhs) ≤ tol
        @test niters(hJAC) == niters(hCG)
    end
end

end


text:
@defcomp damages begin
    DAMFRAC = Variable(index=[time])    # Damages (fraction of gross output)

    TATM    = Parameter(index=[time])   # Increase temperature of atmosphere (degrees C from 1900)
    YGROSS  = Parameter(index=[time])   # Gross world product GROSS of abatement and damages (trillions 2005 USD per year)
    a1      = Parameter()               # Damage intercept
    a2      = Parameter()               # Damage quadratic term
    a3      = Parameter()               # Damage exponent

    TotSLR  = Parameter(index=[time])   # Path of total SLR
    b1      = Parameter()               # Coefficient on SLR
    b2      = Parameter()               # Coefficient on quadratic SLR term
    b3      = Parameter()               # SLR exponent

    function run_timestep(p, v, d, t)
        # Define function for DAMFRAC
        v.DAMFRAC[t] = p.a1 * p.TATM[t] + p.a2 * p.TATM[t]^p.a3 + p.b1 * p.TotSLR[t] + p.b2 * p.TotSLR[t]^p.b3
    end

end


text:
## =========================== ## ## =========================== ##
## 	   Package of collective motion models			 ##
##	   Martin Zumaya Hernandez                               ##
##	   Release Version 0.1.0                                 ##
##     29 / 05 / 2018                                            ##
## =========================== ## ## =========================== ##

module CollectiveDynamics

## SIMPLE VICSEK MODEL 2D
module SVM2D
export Box, Flock, set_output_data_structure_vsk, evolve_system
include("simple_vicsek_model_2D.jl")
end

## SIMPLE VICSEK MODEL 3D
module SVM3D
using Quaternions
export Box, Flock, set_output_data_structure_vsk, evolve_system
include("simple_vicsek_model_3D.jl")
end

## BEHAVIOURAL RULES MODEL
module BehaviouralRules
using Quaternions
export set_output_data_structure, evolve_system
include("behavioural_rules_model.jl")
end

## SHORT AND LONG RANGE INTERACTIONS MODEL
module ShortLongRange
using Distributions, Quaternions
export set_output_data_structure,
    evolve_metric_system,
    evolve_metric_system_2D,
    evolve_topological_system,
    evolve_topological_system_2D
include("short_long_range_model.jl")
end

## INTERIAL SPIN MODEL MODEL
module InertialSpin
using Quaternions
export InertialFlock,
    InertialExtFlock,
    SharedInertialExtFlock,
    InertialParameters,
    set_output_data_structure,
    set_output_data_structure_lr,
    evolve_system,
    evolve_extended_system,
    evolve_extended_system_sh
include("inertial_spin_model.jl")
end

## NETWORK BASED MODEL
module NetworkModel
using Distributions
export set_output_data_structure, set_directed_random_network, evolve_system_2D
include("network_model.jl")
end

end


text:
using SpecialPolynomials;
const SP = SpecialPolynomials;
using Polynomials
import LinearAlgebra: eigvals, det, dot, norm
import SpecialFunctions: gamma
using Test

@testset "Orthogonal" begin
    include("Orthogonal.jl")
end
@testset "Orthogonal compare with tables" begin
    include("Orthogonal-compare.jl")
end
#@testset "ChebyshevT" begin include("ChebyshevT.jl") end
#@testset "ChebyshevU" begin include("ChebyshevU.jl") end
#@testset "Laguerre{α}" begin include("GeneralizedLaguerre.jl") end
#@testset "Gegenbauer{α}" begin include("Gegenbauer.jl") end

@testset "Lagrange" begin
    include("Lagrange.jl")
end
@testset "Newton" begin
    include("Newton.jl")
end

@testset "Bernstein" begin
    include("Bernstein.jl")
end


text:
@testset "I340" begin
    msgs = lintstr("""
    function f(x)
        local a = 1
        local b::Int = 2
        return x+b
    end
    """)
    @test messageset(msgs) == Set([:I340])
    @test contains(msgs[1].message, "unused local variable")
    @test_broken msgs[1].line == 2

    msgs = lintstr("""
    function f(x)
        local (a,b) = (1,2)
        return x+b
    end
    """)
    @test messageset(msgs) == Set([:I340])
    @test contains(msgs[1].message, "unused local variable")
    @test_broken msgs[1].line == 2

    msgs = lintstr("""
    function f(x)
        let a = 1
            b = 2
            y = x+b
            println(y)
        end
    end
    """)
    @test messageset(msgs) == Set([:I340])
    @test contains(msgs[1].message, "unused local variable")

    msgs = lintstr("""
    function f(x)
        local a::Int
        local b = 2::Int # type assertion/conversion
        return x+b
    end
    """)
    @test messageset(msgs) == Set([:I340])
    @test contains(msgs[1].message, "unused local variable")

    msgs = lintstr("""
    function f(x)
        local a
        local b = 2::Int # type assertion/conversion
        return x+b
    end
    """)
    @test messageset(msgs) == Set([:I340])
    @test contains(msgs[1].message, "unused local variable")

    @test isempty(lintstr("""
    function f(x)
        x+=1
    end
    """))

    @test isempty(lintstr("""
    function f(x...)
        Dict(x...)
    end
    """))
end


text:
using Distributed
@everywhere using Mamba

## Data
eyes = Dict{Symbol, Any}(
  :y =>
    [529.0, 530.0, 532.0, 533.1, 533.4, 533.6, 533.7, 534.1, 534.8, 535.3,
     535.4, 535.9, 536.1, 536.3, 536.4, 536.6, 537.0, 537.4, 537.5, 538.3,
     538.5, 538.6, 539.4, 539.6, 540.4, 540.8, 542.0, 542.8, 543.0, 543.5,
     543.8, 543.9, 545.3, 546.2, 548.8, 548.7, 548.9, 549.0, 549.4, 549.9,
     550.6, 551.2, 551.4, 551.5, 551.6, 552.8, 552.9, 553.2],
  :N => 48,
  :alpha => [1, 1]
)

## Model Specification
model = Model(

  y = Stochastic(1,
    (lambda, T, s2, N) ->
      begin
        sigma = sqrt(s2)
        UnivariateDistribution[(
          mu = lambda[Int(T[i])];
          Normal(mu, sigma)) for i in 1:N
        ]
      end,
    false
  ),

  T = Stochastic(1,
    (P, N) -> UnivariateDistribution[Categorical(P) for i in 1:N],
    false
  ),

  P = Stochastic(1,
    alpha -> Dirichlet(alpha)
  ),

  lambda = Logical(1,
    (lambda0, theta) -> Float64[lambda0; lambda0 + theta]
  ),

  lambda0 = Stochastic(
    () -> Normal(0.0, 1000.0),
    false
  ),

  theta = Stochastic(
    () -> Uniform(0.0, 1000.0),
    false
  ),

  s2 = Stochastic(
    () -> InverseGamma(0.001, 0.001)
  )

)

## Initial Values
inits = [
  Dict(:y => eyes[:y], :T => fill(1, eyes[:N]), :P => [0.5, 0.5],
       :lambda0 => 535, :theta => 5, :s2 => 10),
  Dict(:y => eyes[:y], :T => fill(1, eyes[:N]), :P => [0.5, 0.5],
       :lambda0 => 550, :theta => 1, :s2 => 1)
]

## Sampling Scheme
scheme = [DGS(:T),
          Slice([:lambda0, :theta], [5.0, 1.0]),
          Slice(:s2, 2.0, transform=true),
          SliceSimplex(:P, scale=0.75)]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, eyes, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)


text:
"""

    KnetArray{T}(dims)
    KnetArray(a::AbstractArray)
    Array(k::KnetArray)

Container for GPU arrays that supports most of the AbstractArray
interface.  The constructor allocates a KnetArray in the currently
active device, as specified by `gpu()`.  KnetArrays and Arrays can be
converted to each other as shown above, which involves copying to and
from the GPU memory.  Only Float32/64 KnetArrays are fully supported.

Important differences from the alternative CudaArray are: (1) a custom
memory manager that minimizes the number of calls to the slow
cudaMalloc by reusing already allocated but garbage collected GPU
pointers.  (2) a custom getindex that handles ranges such as `a[5:10]`
as views with shared memory instead of copies.  (3) custom CUDA
kernels that implement elementwise, broadcasting, and reduction
operations.

# Supported functions:

* Array operations: ==, !=, cat, convert, copy, copy!, deepcopy,
  display, eachindex, eltype, endof, fill!, first, getindex, hcat,
  isapprox, isempty, length, linearindexing, ndims, ones, pointer,
  rand!, reshape, setindex!, similar, size, stride, strides, summary,
  vcat, vec, zeros.  (Only Integer, Colon, and UnitRange indices
  supported for get/setindex.  CartesianIndex, StepRange, Array, and
  Bool indices not supported.  cat(i,x,y) supported for i=1,2.)

* Math operators: (-), abs, abs2, acos, acosh, asin, asinh, atan,
  atanh, cbrt, ceil, cos, cosh, cospi, erf, erfc, erfcinv, erfcx,
  erfinv, exp, exp10, exp2, expm1, floor, log, log10, log1p, log2,
  round, sign, sin, sinh, sinpi, sqrt, tan, tanh, trunc

* Broadcasting operators: (.*), (.+), (.-), (./), (.<), (.<=), (.!=),
  (.==), (.>), (.>=), (.^), max, min.  (Only Array-Scalar and
  Array-Vector broadcasting are supported. Boolean operators generate
  outputs with same type as inputs; no support for KnetArray{Bool}.)

* Reduction operators: countnz, maximum, minimum, prod, sum, sumabs,
  sumabs2, vecnorm.  (Only Array->Scalar and Array->Vector reductions
  are supported)

* Linear algebra: (*), axpy!, permutedims (only 2D and 3D), transpose

* Knet extras: relu, sigm, invx, logp, logsumexp, conv4, pool,
  deconv4, unpool, mat, update! (Only 4D/5D, Float32/64 KnetArrays
  support conv4, pool, deconv4, unpool)

# Memory management

Knet models do not overwrite arrays which need to be preserved for
gradient calculation.  This leads to a lot of allocation and regular
GPU memory allocation is prohibitively slow. Fortunately most models
use identically sized arrays over and over again, so we can minimize
the number of actual allocations by reusing preallocated but garbage
collected pointers.

When Julia gc reclaims a KnetArray, a special finalizer keeps its
pointer in a table instead of releasing the memory.  If an array with
the same size in bytes is later requested, the same pointer is reused.
The exact algorithm for allocation is:

1. Try to find a previously allocated and garbage collected pointer in
   the current device. (0.5 μs)

2. If not available, try to allocate a new array using cudaMalloc. (10
   μs)

3. If not successful, try running gc() and see if we get a pointer of
   the right size. (75 ms, but this should be amortized over all
   reusable pointers that become available due to the gc)

4. Finally if all else fails, clean up all saved pointers in the
   current device using cudaFree and try allocation one last
   time. (25-70 ms, however this causes the elimination of all
   reusable pointers)

"""
type KnetArray{T,N}
    ptr::KnetPtr
    dims::NTuple{N,Int}
end

# Note: I removed <: AbstractArray{T,N} after I painfully discovered
# some inefficient AbstractArray methods inherited unintentionally.
# It is better to define a few extra methods to keep a tighter control
# on what methods exactly get called for KnetArrays.

# Aliases:
typealias KnetMatrix{T} KnetArray{T,2}
typealias KnetVector{T} KnetArray{T,1}
typealias KnetVecOrMat{T} Union{KnetVector{T}, KnetMatrix{T}}

# Constructors:
import Base: convert
# Internal constructor defines KnetArray{T,N}(ptr,dims)
# These define KnetArray{T,N}(dims) and KnetArray{T,N}(d...)
if VERSION >= v"0.5.0-dev+7720"; @eval begin
    (::Type{KnetArray{T,N}}){T,N}(d::NTuple{N,Int})=KnetArray{T,N}(KnetPtr(sizeof(T)*prod(d)), d)
    (::Type{KnetArray{T,N}}){T,N}(d::Int...) = KnetArray{T,N}(d)
    (::Type{KnetArray{T,N}}){T,N}(d::Integer...) = KnetArray{T,N}(convert(Tuple{Vararg{Int}}, d))
end; else; @eval begin
    convert{T,N}(::Type{KnetArray{T,N}},d::NTuple{N,Int})=KnetArray{T,N}(KnetPtr(sizeof(T)*prod(d)), d)
    convert{T,N}(::Type{KnetArray{T,N}},d::Int...)=convert(KnetArray{T,N},d)
    convert{T,N}(::Type{KnetArray{T,N}},d::Integer...)=convert(KnetArray{T,N},convert(Tuple{Vararg{Int}}, d))
end; end
# These define KnetArray{T}(dims) and KnetArray{T}(d...)
if VERSION >= v"0.5.0-dev+7720"; @eval begin
    (::Type{KnetArray{T}}){T,N}(d::NTuple{N,Int}) = KnetArray{T,N}(d)
    (::Type{KnetArray{T}}){T}(d::Int...) = KnetArray{T}(d)
    (::Type{KnetArray{T}}){T}(d::Integer...) = KnetArray{T}(convert(Tuple{Vararg{Int}}, d))
end; else; @eval begin
    convert{T,N}(::Type{KnetArray{T}},d::NTuple{N,Int})=KnetArray{T,N}(KnetPtr(sizeof(T)*prod(d)), d)
    convert{T}(::Type{KnetArray{T}},d::Int...)=KnetArray{T}(d)
    convert{T}(::Type{KnetArray{T}},d::Integer...)=KnetArray{T}(convert(Tuple{Vararg{Int}}, d))
end; end
# These define KnetArray(T,dims) and KnetArray(T,d...)
KnetArray{T,N}(::Type{T}, d::NTuple{N,Int}) = KnetArray{T,N}(d)
KnetArray{T}(::Type{T}, d::Int...)=KnetArray(T,d)
KnetArray{T}(::Type{T}, d::Integer...)=KnetArray(T,convert(Tuple{Vararg{Int}},d))

# Conversions:
import Base: convert, reshape, vec, unsafe_convert, pointer
# KnetArray <- KnetArray
convert{T,N}(::Type{KnetArray}, x::KnetArray{T,N}) = x
convert{T,N}(::Type{KnetArray{T}}, x::KnetArray{T,N}) = x
convert{T,N}(::Type{KnetArray{T,N}}, x::KnetArray{T,N}) = x
convert{T,N,S}(::Type{KnetArray{T}}, x::KnetArray{S,N}) = convert(KnetArray{T,N}, x)
convert{T,N,S}(::Type{KnetArray{T,N}}, x::KnetArray{S,N}) = convert(KnetArray{T,N},unsafe_copy!(Array(S, size(x)), 1, x, 1, length(x)))
reshape{T}(a::KnetArray{T},dims::Dims)=(if dims==size(a); a; elseif prod(dims)!=length(a); throw(DimensionMismatch()); else; KnetArray{T,length(dims)}(a.ptr,dims); end)
reshape(a::KnetArray, dims::Int...) = reshape(a, dims)
vec(a::KnetArray) = reshape(a, length(a))
# KnetArray <- AbstractArray
convert{T,N}(::Type{KnetArray}, x::AbstractArray{T,N}) = convert(KnetArray{T,N}, x)
convert{T,N,S}(::Type{KnetArray{T}}, x::AbstractArray{S,N}) = convert(KnetArray{T,N}, x)
convert{T,N,S}(::Type{KnetArray{T,N}}, x::AbstractArray{S,N}) = unsafe_copy!(KnetArray{T}(size(x)), 1, convert(Array{T,N},x), 1, length(x))
# Array <- KnetArray
convert{T,N}(::Type{Array}, x::KnetArray{T,N}) = convert(Array{T,N}, x)
convert{T,N,S}(::Type{Array{T}}, x::KnetArray{S,N}) = convert(Array{T,N}, x)
convert{T,N,S}(::Type{Array{T,N}}, x::KnetArray{S,N}) = convert(Array{T,N},unsafe_copy!(Array{S}(size(x)), 1, x, 1, length(x)))
# Ptr <- KnetArray
unsafe_convert{T}(::Type{Ptr{T}}, a::KnetArray) = unsafe_convert(Ptr{T}, pointer(a))
pointer{T}(a::KnetArray{T})=convert(Ptr{T}, a.ptr.ptr)
pointer{T}(a::KnetArray{T},i)=convert(Ptr{T}, a.ptr.ptr + (i-1)*sizeof(T))

# AbstractArray interface
import Base: eachindex, eltype, endof, fill!, first, isempty, length, linearindexing, ndims, ones, similar, size, stride, strides, zeros, (==), isapprox
eachindex(a::KnetArray) = (1:length(a))
eltype{T}(::KnetArray{T})=T
eltype{T}(::Type{KnetArray{T}}) = T
eltype{T,n}(::Type{KnetArray{T,n}}) = T
endof(a::KnetArray) = length(a)
fill!{T}(a::KnetArray{T},x)=(knetfill!(a,T(x),1,length(a));a)
first(a::KnetArray) = a[1]
# AutoGrad leaves `first` as a compound proc calling start which doesn't work with KnetArrays
@primitive  first(x::KnetArray),dy,y  AutoGrad.ungetindex(dy,x,1)
isempty(a::KnetArray) = (0==length(a))
length(a::KnetArray)=prod(size(a))
linearindexing(::KnetArray)=Base.LinearFast()
ndims{T,N}(a::KnetArray{T,N})=N
ones{T}(a::KnetArray{T})=fill!(similar(a),one(T))
similar(a::KnetArray, T, dims::Dims)      = KnetArray(T, dims)
similar(a::KnetArray, T, dims::Int...)    = similar(a, T, dims)
similar(a::KnetArray, T)                  = similar(a, T, size(a))
similar{T}(a::KnetArray{T})               = similar(a, T, size(a))
similar{T}(a::KnetArray{T}, dims::Dims)   = similar(a, T, dims)
similar{T}(a::KnetArray{T}, dims::Int...) = similar(a, T, dims)
size(a::KnetArray)=a.dims
size{T,N}(a::KnetArray{T,N},i::Integer)=(if i>N; 1; else; size(a)[i]; end)
stride{T,N}(a::KnetArray{T,N},i::Integer)=(if i>N; length(a); else; s=1; for n=1:(i-1); s*=size(a,n); end; s; end)
strides{T,N}(a::KnetArray{T,N})=ntuple(n->stride(a,n), N)
zeros{T}(a::KnetArray{T})=fill!(similar(a),zero(T))

# Comparisons
(==){T}(a::KnetArray{T},b::KnetArray{T})=(size(a)==size(b) && vecnorm(a-b)==0)
(==)(a::AbstractArray,b::KnetArray)=(size(a)==size(b) && a==Array(b))
(==)(a::KnetArray,b::AbstractArray)=(size(a)==size(b) && Array(a)==b)
# Adapted from base/linalg/generic.jl:589
isapprox{T}(a::KnetArray{T}, b::KnetArray{T}; rtol=sqrt(eps(T)), atol=T(0))=(size(a)==size(b) && vecnorm(a-b) <= atol + rtol * max(vecnorm(a), vecnorm(b)))
isapprox(a::AbstractArray,b::KnetArray;o...)=(size(a)==size(b) && isapprox(a,Array(b);o...))
isapprox(a::KnetArray,b::AbstractArray;o...)=(size(a)==size(b) && isapprox(Array(a),b;o...))

# Indexing:
import Base: getindex, setindex!

# We will implement indexing ranges as views not copies, if possible (when contiguous).
# For contiguous memory without stride all but the last >1 dimension must be full

# The original getindex(a,i:j...) for AbstractArray copies:
# function _getindex(l::LinearIndexing, A::AbstractArray, I::Union{Real, AbstractArray, Colon}...)
# in abstractarray.jl:487,multidimensional.jl:184.

# which getindex ops does array implement?
# getindex(A::Array, i1::Real)
# getindex(A::Array, i1::Real, i2::Real, I::Real...)
# getindex(A::Array, I::UnitRange{Int})
# getindex(A::Array, c::Colon)
# getindex{T<:Real}(A::Array, I::Range{T})

# Julia #14770
# If I is shorter than ndims(A) but longer than 1 the remaining indices assumed =1
# Also extra 1's at the end of I are ignored

# These two are not sufficient in spite of what the documentation says:
# display goes into an infinite loop!
# getindex{T}(A::KnetArray{T}, i::Int)=unsafe_copy!(T[0], 1, A, i, 1)[1]
# setindex!{T}(A::KnetArray{T}, v, i::Int)=unsafe_copy!(A, i, T[v], 1, 1)

# First deal with the easy cases: integer indices, a Colon or a UnitRange.

function getindex{T}(A::KnetArray{T}, I::Real)
    J = Int(I)
    1 <= J <= length(A) || throw(BoundsError(A,J))
    unsafe_copy!(T[0], 1, A, J, 1)[1]
end

function setindex!{T}(A::KnetArray{T}, v, I::Real)
    J = Int(I)
    1 <= J <= length(A) || throw(BoundsError(A,J))
    unsafe_copy!(A, J, T[v], 1, 1)
end

function getindex{T}(A::KnetArray{T}, I::Real...)
    J = Base.to_indexes(I...)
    @inbounds for j=1:length(J)
        1 <= J[j] <= size(A,j) || throw(BoundsError(A,J))
    end
    i = sub2ind(size(A), J...)
    unsafe_copy!(T[0], 1, A, i, 1)[1]
end

function setindex!{T}(A::KnetArray{T}, v, I::Real...)
    J = Base.to_indexes(I...)
    @inbounds for j=1:length(J)
        1 <= J[j] <= size(A,j) || throw(BoundsError(A,J))
    end
    i = sub2ind(size(A), J...)
    unsafe_copy!(A, i, T[v], 1, 1)
end

function getindex{T}(A::KnetArray{T}, I::UnitRange)
    1 <= first(I) <= last(I) <= length(A) || throw(BoundsError(A,I))
    off = 1+(first(I)-1)*sizeof(T)
    len = length(I)*sizeof(T)
    ptr = KnetPtr(A.ptr, off, len)
    KnetArray{T,1}(ptr, (length(I),))
end

function setindex!{T}(A::KnetArray{T}, v, I::UnitRange)
    1 <= first(I) <= last(I) <= length(A) || throw(BoundsError(A,I))
    if isa(v,Number)
        knetfill!(A,T(v),first(I),length(I))
    elseif (isa(v,KnetArray) || isa(v,Array))
        length(v)==length(I) || throw(DimensionMismatch())
        eltype(v)==T || (v = convert(Array{T},v))
        unsafe_copy!(A,first(I),v,1,length(I))
    else
        throw(MethodError(setindex!, A, v, I))
    end
end

function getindex(A::KnetArray, I::Colon)
    reshape(A,length(A))
end

function setindex!{T}(A::KnetArray{T}, v, I::Colon)
    if isa(v,Number)
        knetfill!(A,T(v),1,length(A))
    elseif (isa(v,KnetArray) || isa(v,Array))
        length(v)==length(A) || throw(DimensionMismatch())
        eltype(v)==T || (v = convert(Array{T},v))
        unsafe_copy!(A,1,v,1,length(A))
    else
        throw(MethodError(setindex!, A, v, I))
    end
end

# TODO: the following getindex, setindex! work for 1 and 2 dimensions only, write general versions.

function getindex{T,N}(A::KnetArray{T,N}, I::Union{Real, UnitRange, Colon}...)
    (nelts,nrows,ncols,firstindex,astep) = indexparams(A,I...)
    B1 = isa(I[1],Colon) ? size(A,1) : length(I[1])
    B2 = isa(I[2],Colon) ? size(A,2) : length(I[2])
    if ncols == 1
        off = 1+(firstindex-1)*sizeof(T)
        len = nrows*sizeof(T)
        ptr = KnetPtr(A.ptr, off, len)
        KnetArray{T,2}(ptr, (B1,B2))
    else
        B = similar(A, (B1,B2))
        nrows *= sizeof(T); astep *= sizeof(T)
        ccall((:xcopy,libknet8),Void,(Cint,Cint,Cptr,Cint,Cptr,Cint),
              nrows, ncols, pointer(A,firstindex), astep, B, nrows)
        return B
    end
end

function setindex!{T,N}(A::KnetArray{T,N}, B, I::Union{Real, UnitRange, Colon}...)
    (nelts,nrows,ncols,firstindex,astep) = indexparams(A,I...)
    aptr0 = pointer(A, firstindex)
    if isa(B,Number)
        B = T(B)
        if ncols == 1
            if T <: Float32;    ccall((:fill_32,libknet8),Void,(Cint,Cfloat, Ptr{Cfloat}), nelts,B,aptr0)
            elseif T<: Float64; ccall((:fill_64,libknet8),Void,(Cint,Cdouble,Ptr{Cdouble}),nelts,B,aptr0)
            else error("$T not supported"); end
        else
            if T <: Float32;    ccall((:xfill_32,libknet8),Void,(Cint,Cint,Cfloat, Ptr{Cfloat}, Cint),nrows,ncols,B,aptr0,astep)
            elseif T<: Float64; ccall((:xfill_64,libknet8),Void,(Cint,Cint,Cdouble,Ptr{Cdouble},Cint),nrows,ncols,B,aptr0,astep)
            else error("$T not supported"); end
        end
    else
        length(B) == nelts || throw(DimensionMismatch())
        B = convert(KnetArray{T},B)
        if ncols == 1
            @cuda(cudart,cudaMemcpyAsync,(Cptr,Cptr,Csize_t,UInt32,Cptr),
                  aptr0, B, nelts*sizeof(T), cudadir(A,B), C_NULL)
        else
            nrows *= sizeof(T); astep *= sizeof(T)
            ccall((:xcopy,libknet8),Void,(Cint,Cint,Cptr,Cint,Cptr,Cint), nrows, ncols, B, nrows, aptr0, astep)
        end
    end
    return A
end

function indexparams{T,N}(A::KnetArray{T,N}, I::Union{Real, UnitRange, Colon}...)
    N > 2 && error("setindex for ndims > 2 not implemented yet")
    skipped = false
    nrows = nelts = 1
    subs1 = ones(Int,2)
    astep = length(A)
    for i=1:length(I)
        Ii = I[i]
        Ai = size(A,i)
        if isa(Ii, Colon)
            Li = Ai
            subs1[i] = 1
        elseif isa(Ii, Real)
            1 <= Ii <= Ai || throw(DimensionMismatch())
            Li = 1
            subs1[i] = Int(Ii)
        else
            1 <= first(Ii) <= last(Ii) <= Ai || throw(DimensionMismatch())
            Li = length(Ii)
            subs1[i] = first(Ii)
        end
        nelts *= Li
        if !skipped
            nrows *= Li
            if Li < Ai
                skipped = true
                astep = stride(A,i+1)
            end
        end
    end
    ncols = div(nelts, nrows)
    firstindex = sub2ind(size(A),subs1...)
    return (nelts,nrows,ncols,firstindex,astep)
end

# Concatenation:
import Base: hcat, vcat, cat

# Benchmarks in μs for hcat and vcat: a=rand(1000,1000) v=rand(1000), t=v'
#		cpu	gpu	g->c->g	vkernel
# hcat(a,a)	2350	225	16160
# hcat(a,v)	1230	115	6490
# hcat(v,a)	1220	120	6490
# hcat(v,v)	3.53	12.53	48.49
# vcat(a,a)	2630	10980	16590	665
# vcat(a,t)	1350	10860	6550	338
# vcat(t,a)	1360	10850	6570	338
# vcat(v,v)	2.13	12.33	45.40	13.58

# setindex! methods called by hcat/vcat:
# hcat(v,v): I = (Colon(),1:1) I = (Colon(),2:2)
# vcat(v,v): uses single index
# hcat(m,m): I = (Colon(),1:5) I = (Colon(),6:10)
# vcat(m,m): I = (1:3,Colon()) I = (4:6,Colon())

function hcat{T}(a::KnetVecOrMat{T}, b::KnetVecOrMat{T})
    size(a,1)==size(b,1) || throw(DimensionMismatch())
    c1 = size(a,1)
    c2 = size(a,2) + size(b,2)
    c = KnetArray(T, (c1,c2))
    c[:,1:size(a,2)] = a
    c[:,1+size(a,2):end] = b
    return c
end

function vcat{T}(a::KnetVector{T}, b::KnetVector{T})
    c = KnetArray(T, length(a)+length(b))
    c[1:length(a)] = a
    c[1+length(a):end] = b
    return c
end

function vcat{T}(a::KnetVecOrMat{T}, b::KnetVecOrMat{T})
    size(a,2)==size(b,2) || throw(DimensionMismatch())
    c1 = size(a,1) + size(b,1)
    c2 = size(a,2)
    c = KnetArray(T, (c1,c2))
    c[1:size(a,1),:] = a
    c[1+size(a,1):end,:] = b
    return c
end

function cat{T}(d, a::KnetVecOrMat{T}, b::KnetVecOrMat{T})
    if     d==1; vcat(a,b)
    elseif d==2; hcat(a,b)
    else error("cat($d) not implemented.")
    end
end

# Utilities:

# Generalizing low level copy using linear indexing to/from gpu arrays:
# copy!{T}(dest::Array{T}, doffs::Integer, src::Array{T}, soffs::Integer, n::Integer)
# Note that this is an unsafe operation, no argument or bounds checking performed.
# Defined in Base:
# unsafe_copy!{T}(dest::Ptr{T}, src::Ptr{T}, n) at array.jl:73
# unsafe_copy!{T}(dest::Array{T,N}, doffs, src::Array{T,N}, soffs, n) at array.jl:79

import Base: unsafe_copy!, copy, copy!
typealias KorA{T} Union{KnetArray{T},Array{T}}

function copy!{T}(dest::KorA{T}, doffs::Integer, src::KorA{T}, soffs::Integer, n::Integer; stream=C_NULL)
    n == 0 && return dest
    n > 0 || throw(ArgumentError(string("tried to copy n=", n, " elements, but n should be nonnegative")))
    if soffs < 1 || doffs < 1 || soffs+n-1 > length(src) || doffs+n-1 > length(dest)
        throw(BoundsError())
    end
    unsafe_copy!(dest, doffs, src, soffs, n; stream=stream)
end

copy!{T}(dest::KorA{T}, src::KorA{T}) = copy!(dest, 1, src, 1, length(src))

copy(a::KnetArray)=unsafe_copy!(similar(a),1,a,1,length(a))

# This will make deepcopy work properly
Base.deepcopy_internal(x::KnetArray, s::ObjectIdDict)=if haskey(s,x); s[x]; else; copy(x); end

function unsafe_copy!{T}(dest::KorA{T}, doffs, src::KorA{T}, soffs, n; stream=C_NULL)
    @cuda(cudart,cudaMemcpyAsync,(Cptr,Cptr,Csize_t,UInt32,Cptr),
          pointer(dest,doffs), pointer(src,soffs), n*sizeof(T), cudadir(dest,src), stream)
    return dest
end

function cudadir(a,b)
    deva = isa(a,KnetArray) && a.ptr.dev >= 0
    devb = isa(b,KnetArray) && b.ptr.dev >= 0
    if !deva && !devb; return 0
    elseif deva && !devb; return 1
    elseif !deva && devb; return 2
    elseif deva && devb;  return 3
    end
end

# Efficient fill:
for S in (32,64); T = Symbol("Float$S"); F = "fill_$S"
    @eval function knetfill!(a::KnetArray{$T},v::$T,off,len)
        ccall(($F,$libknet8),Void,(Cint,$T,Ptr{$T}),len,v,pointer(a,off))
    end
end

# AutoGrad functions:
import AutoGrad: zeroslike, sum_outgrads, OneHot, unary_nd, indexed_function, isequivalent
zeroslike(a::KnetArray)=zeros(a)
sum_outgrads{T}(a::KnetArray{T},b::KnetArray{T})=(a+b)
sum_outgrads(a::KnetArray,b::OneHot)=setindex!(a,sum_outgrads(getindex(a,b.index...),b.value),b.index...)
unary_nd(f, x::KnetArray, eps) = reshape(eltype(x)[unary_nd(indexed_function(f, x, i), x[i], eps) for i in 1:length(x)], size(x))
isequivalent(x::Union{KnetArray,AbstractArray}, y::Union{KnetArray,AbstractArray}; o...)=(length(x)==length(y) && all(i->isequivalent(x[i],y[i];o...), 1:length(x)))

# Hack for printing without copying the whole KnetArray and without inheriting AbstractArray:
import Base: display, summary
type KnetDisplay{T,N} <: AbstractArray{T,N}; a::KnetArray{T,N}; end
getindex(a::KnetDisplay, i...) = getindex(a.a, i...)
size(a::KnetDisplay) = size(a.a)
summary(a::KnetDisplay) = summary(a.a)
summary(a::KnetArray) = string(Base.dims2string(size(a)), " ", typeof(a))
display(a::KnetArray) = display(KnetDisplay(a))
AutoGrad._dbg(a::KnetArray) = "K$(join([AutoGrad.id2(a),size(a)...],'_'))"

# curand functions:

import Base: rand!
rand!(a::KnetArray{Float32})=(@cuda(curand,curandGenerateUniform,(Cptr,Ptr{Cfloat},Csize_t),rng(),a,length(a)); a)
rand!(a::KnetArray{Float64})=(@cuda(curand,curandGenerateUniformDouble,(Cptr,Ptr{Cdouble},Csize_t),rng(),a,length(a)); a)

let RNG=0
global rng
function rng(init=false)
    if RNG==0 || init
        ptr = Cptr[0]
        # CURAND_RNG_PSEUDO_DEFAULT = 100, ///< Default pseudorandom generator
        @cuda(curand,curandCreateGenerator,(Cptr,Cint),ptr,100)
        RNG = ptr[1]
    end
    return RNG
end
end

# Array/KnetArray Transfer

# This works but unnecessarily defines new functions:
# cpu2gpu(x::Array)=KnetArray(x)
# @primitive cpu2gpu(x),dy,y (gpu2cpu(dy))
# gpu2cpu(x::KnetArray)=Array(x)
# @primitive gpu2cpu(x),dy,y (cpu2gpu(dy))

# This does not work because !isa(Array,Function)
# @primitive  KnetArray(x::Array),dy  Array(dy)
# @primitive  Array(x::KnetArray),dy  KnetArray(dy)

# This does not work, parametric methods not yet supported, also unnecessary first arg gradient.
# @primitive convert{A<:AbstractArray,K<:KnetArray}(T::Type{K}, x::Rec{A}),dy 0 Array(dy)
# @primitive convert{A<:AbstractArray,K<:KnetArray}(T::Type{A}, x::Rec{K}),dy 0 KnetArray(dy)

# So we will define gradients for convert, KnetArray, Array manually:
Base.Array{K<:KnetArray}(x::Rec{K})=convert(Array,x)
KnetArray{A<:AbstractArray}(x::Rec{A})=convert(KnetArray,x)
let convert_r = recorder(convert)
    global convert
    convert(::Type{Grad{2}},dy,y,T,x) = convert(typeof(AutoGrad.getval(x)),dy)
    # This does not work, it breaks the Node(::Rec) constructor, so we define Knet specific version.
    # convert(T::Type, x::Rec) = convert_r(T,x)
    convert{A<:AbstractArray,K<:KnetArray}(::Type{A},x::Rec{K})=convert_r(A,x)
    convert{A<:AbstractArray,K<:KnetArray}(::Type{K},x::Rec{A})=convert_r(K,x)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1VolumeMount <: SwaggerModel
    mountPath::Any # spec type: Union{ Nothing, String } # spec name: mountPath
    mountPropagation::Any # spec type: Union{ Nothing, String } # spec name: mountPropagation
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    subPath::Any # spec type: Union{ Nothing, String } # spec name: subPath
    subPathExpr::Any # spec type: Union{ Nothing, String } # spec name: subPathExpr

    function IoK8sApiCoreV1VolumeMount(;mountPath=nothing, mountPropagation=nothing, name=nothing, readOnly=nothing, subPath=nothing, subPathExpr=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1VolumeMount, Symbol("mountPath"), mountPath)
        setfield!(o, Symbol("mountPath"), mountPath)
        validate_property(IoK8sApiCoreV1VolumeMount, Symbol("mountPropagation"), mountPropagation)
        setfield!(o, Symbol("mountPropagation"), mountPropagation)
        validate_property(IoK8sApiCoreV1VolumeMount, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiCoreV1VolumeMount, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1VolumeMount, Symbol("subPath"), subPath)
        setfield!(o, Symbol("subPath"), subPath)
        validate_property(IoK8sApiCoreV1VolumeMount, Symbol("subPathExpr"), subPathExpr)
        setfield!(o, Symbol("subPathExpr"), subPathExpr)
        o
    end
end # type IoK8sApiCoreV1VolumeMount

const _property_map_IoK8sApiCoreV1VolumeMount = Dict{Symbol,Symbol}(Symbol("mountPath")=>Symbol("mountPath"), Symbol("mountPropagation")=>Symbol("mountPropagation"), Symbol("name")=>Symbol("name"), Symbol("readOnly")=>Symbol("readOnly"), Symbol("subPath")=>Symbol("subPath"), Symbol("subPathExpr")=>Symbol("subPathExpr"))
const _property_types_IoK8sApiCoreV1VolumeMount = Dict{Symbol,String}(Symbol("mountPath")=>"String", Symbol("mountPropagation")=>"String", Symbol("name")=>"String", Symbol("readOnly")=>"Bool", Symbol("subPath")=>"String", Symbol("subPathExpr")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1VolumeMount }) = collect(keys(_property_map_IoK8sApiCoreV1VolumeMount))
Swagger.property_type(::Type{ IoK8sApiCoreV1VolumeMount }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1VolumeMount[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1VolumeMount }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1VolumeMount[property_name]

function check_required(o::IoK8sApiCoreV1VolumeMount)
    (getproperty(o, Symbol("mountPath")) === nothing) && (return false)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1VolumeMount }, name::Symbol, val)
end


text:
module MyApp

using Example
using Pkg.Artifacts

socrates = joinpath(ensure_artifact_installed("socrates", joinpath(@__DIR__, "..", "Artifacts.toml")),
                    "bin", "socrates")

Base.@ccallable function julia_main()::Cint
    try
        real_main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

function real_main()
    @show ARGS
    @show Base.PROGRAM_FILE
    @show DEPOT_PATH
    @show LOAD_PATH
    @show pwd()
    @show Base.active_project()
    @show Threads.nthreads()
    @show Sys.BINDIR
    display(Base.loaded_modules)
    println("Running an artifact:")
    run(`$socrates`)
    println()
    @show unsafe_string(Base.JLOptions().image_file)
    @show Example.domath(5)
    @show sin(0.0)
    return
end

if abspath(PROGRAM_FILE) == @__FILE__
    real_main()
end

end # module


text:
module TrickyArithmetic
    struct A
        x::Int
    end
    A(a::A) = a
    Base.convert(::Type{A}, i::Int) = A(i)
    Base.zero(::Union{A, Type{A}}) = A(0)
    Base.one(::Union{A, Type{A}}) = A(1)
    struct B
        x::Int
    end
    struct C
        x::Int
    end
    C(a::A) = C(a.x)
    Base.zero(::Union{C, Type{C}}) = C(0)
    Base.one(::Union{C, Type{C}}) = C(1)

    Base.:(*)(x::Int, a::A) = B(x*a.x)
    Base.:(*)(a::A, x::Int) = B(a.x*x)
    Base.:(*)(a::Union{A,B}, b::Union{A,B}) = B(a.x*b.x)
    Base.:(*)(a::Union{A,B,C}, b::Union{A,B,C}) = C(a.x*b.x)
    Base.:(+)(a::Union{A,B,C}, b::Union{A,B,C}) = C(a.x+b.x)
    Base.:(-)(a::Union{A,B,C}, b::Union{A,B,C}) = C(a.x-b.x)

    struct D{NT, DT}
        n::NT
        d::DT
    end
    D{NT, DT}(d::D{NT, DT}) where {NT, DT} = d # called by oneunit
    Base.zero(::Union{D{NT, DT}, Type{D{NT, DT}}}) where {NT, DT} = zero(NT) / one(DT)
    Base.one(::Union{D{NT, DT}, Type{D{NT, DT}}}) where {NT, DT} = one(NT) / one(DT)
    Base.convert(::Type{D{NT, DT}}, a::Union{A, B, C}) where {NT, DT} = NT(a) / one(DT)
    #Base.convert(::Type{D{NT, DT}}, a::D) where {NT, DT} = NT(a.n) / DT(a.d)

    Base.:(*)(a::D, b::D) = (a.n*b.n) / (a.d*b.d)
    Base.:(*)(a::D, b::Union{A,B,C}) = (a.n * b) / a.d
    Base.:(*)(a::Union{A,B,C}, b::D) = b * a
    Base.inv(a::Union{A,B,C}) = A(1) / a
    Base.inv(a::D) = a.d / a.n
    Base.:(/)(a::Union{A,B,C}, b::Union{A,B,C}) = D(a, b)
    Base.:(/)(a::D, b::Union{A,B,C}) = a.n / (a.d*b)
    Base.:(/)(a::Union{A,B,C,D}, b::D) = a * inv(b)
    Base.:(+)(a::Union{A,B,C}, b::D) = (a*b.d+b.n) / b.d
    Base.:(+)(a::D, b::Union{A,B,C}) = b + a
    Base.:(+)(a::D, b::D) = (a.n*b.d+a.d*b.n) / (a.d*b.d)
    Base.:(-)(a::Union{A,B,C}) = typeof(a)(a.x)
    Base.:(-)(a::D) = (-a.n) / a.d
    Base.:(-)(a::Union{A,B,C,D}, b::Union{A,B,C,D}) = a + (-b)

    Base.promote_rule(::Type{A}, ::Type{B}) = B
    Base.promote_rule(::Type{B}, ::Type{A}) = B
    Base.promote_rule(::Type{A}, ::Type{C}) = C
    Base.promote_rule(::Type{C}, ::Type{A}) = C
    Base.promote_rule(::Type{B}, ::Type{C}) = C
    Base.promote_rule(::Type{C}, ::Type{B}) = C
    Base.promote_rule(::Type{D{NT,DT}}, T::Type{<:Union{A,B,C}}) where {NT,DT} = D{promote_type(NT,T),DT}
    Base.promote_rule(T::Type{<:Union{A,B,C}}, ::Type{D{NT,DT}}) where {NT,DT} = D{promote_type(NT,T),DT}
    Base.promote_rule(::Type{D{NS,DS}}, ::Type{D{NT,DT}}) where {NS,DS,NT,DT} = D{promote_type(NS,NT),promote_type(DS,DT)}
end

max_stars_repo_name: perrutquist/julia

max_issues_repo_name: perrutquist/julia

max_forks_repo_name: perrutquist/julia


text:
                 ##################################
    ###################     Muffin Demo     #################
                 ##################################

################################

myfolder = string(Pkg.dir("Muffin"))
mydataobj =  "data/M31.fits"
mydatapsf =  "data/meerkat_m30_25pix.psf.fits"
mynitermax = 5
myrhop = 1
myrhot = 5
myrhov = 2
myrhos = 1
myμt = 5e-1
myμv = 1
mymueps = 1e-3
mysavepath = string(myfolder,myfolder[1],"data/demo_results.jld")

################################

psfst, skyst, algost, admmst, toolst = muffin(folder = myfolder,
                                              dataobj = mydataobj,
                                              datapsf = mydatapsf,
                                              nitermax = mynitermax,
                                              rhop = myrhop,
                                              rhot = myrhot,
                                              rhov = myrhov,
                                              rhos = myrhos,
                                              μt = myμt,
                                              μv = myμv,
                                              mueps = mymueps)

################################

savedata(mysavepath, psfst, skyst, algost, admmst, toolst)

################################


text:
export buffer_remove_row

"""Shifts Array Up rather then allocating a new array
   Returns the new array length (DO NOT USE length(array) after using this function as it is useless)

   TODO: Generialized Solution
"""
function buffer_remove_row(array::AbstractArray{T, N}, removable_parts, true_array_len) where T where N
      if length(removable_parts) == 0; return (array, true_array_len); end
      dimensions = N
      placement_index = minimum(removable_parts)
      index_vector = collect(placement_index:true_array_len)
      for r in 1:length(removable_parts)
            rem = removable_parts[r]
            if rem <= true_array_len
                  deleteat!(index_vector, rem - placement_index - r + 2)
            end
      end
      new_len = size(array, 1) - length(removable_parts)
      for i in index_vector
            if N == 1
                  array[placement_index] = array[i]
            elseif N == 2
                  array[placement_index, 1:end] = array[i, 1:end]
            elseif N == 3
                  array[placement_index, 1:end, 1:end] = array[i, 1:end, 1:end]
            elseif N == 4
                  array[placement_index, 1:end, 1:end, 1:end] = array[i, 1:end, 1:end, 1:end]
            end
            placement_index += 1
      end
      return (array, new_len)
end


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	minimum_capacity_requirement(EP::Model, inputs::Dict)
The minimum capacity requirement constraint allows for modeling minimum deployment of a certain technology or set of eligible technologies across the eligible model zones and can be used to mimic policies supporting specific technology build out (i.e. capacity deployment targets/mandates for storage, offshore wind, solar etc.). The default unit of the constraint is in MW. For each requirement $p \in \mathcal{P}^{MinCapReq}$, we model the policy with the following constraint.
```math
\begin{aligned}
\sum_{y \in \mathcal{G} } \sum_{z \in \mathcal{Z} } \left( \epsilon_{y,z,p}^{MinCapReq} \times \Delta^{\text{total}}_{y,z} \right) \geq REQ_{p}^{MinCapReq} \hspace{1 cm}  \forall p \in \mathcal{P}^{MinCapReq}
\end{aligned}
```
Note that $\epsilon_{y,z,p}^{MinCapReq}$ is the eligiblity of a generator of technology $y$ in zone $z$ of requirement $p$ and will be equal to $1$ for eligible generators and will be zero for ineligible resources. The dual value of each minimum capacity constraint can be interpreted as the required payment (e.g. subsidy) per MW per year required to ensure adequate revenue for the qualifying resources.
"""
function minimum_capacity_requirement(EP::Model, inputs::Dict, setup::Dict)

	println("Minimum Capacity Requirement Module")

	NumberOfMinCapReqs = inputs["NumberOfMinCapReqs"]

	@constraint(EP, cZoneMinCapReq[mincap = 1:NumberOfMinCapReqs], EP[:eMinCapRes][mincap] >= inputs["MinCapReq"][mincap])

	return EP
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
using DBInterface

@info "Loading SQLite support for Material items into NeXLCore."

"""
   buildMaterialTables(db::SQLite.DB)

Build the necessary tables in a `SQLite` database to hold `Material`s.
Safe! Only creates the tables if they don't exist.
"""
function buildMaterialTables(db::SQLite.DB)
    commands = (
        """
        CREATE TABLE IF NOT EXISTS MATERIAL (
            PKEY INTEGER PRIMARY KEY AUTOINCREMENT,
            NAME TEXT NOT NULL, -- Make material name unique
            DENSITY REAL
        );""",
        "CREATE INDEX IF NOT EXISTS MATNAME_INDEX ON MATERIAL(NAME);",
        """
        CREATE TABLE IF NOT EXISTS MASSFRACTION (
            PKEY INTEGER PRIMARY KEY AUTOINCREMENT,
            MATKEY INTEGER NOT NULL,
            MFZ INTEGER NOT NULL, -- Atomic number
            MFC REAL NOT NULL, -- Mass fraction
            MFUC REAL, -- Uncertainty in c (1 σ) or NULL if unknown
            MFA REAL, -- atomic weight or NULL if default
            FOREIGN KEY(MATKEY) REFERENCES MATERIAL(PKEY)
        );""",
        """
        CREATE TABLE IF NOT EXISTS MATERIALPROPERTY ( -- Holds text properties
            MATKEY INTEGER NOT NULL,
            KEY TEXT NOT NULL,
            VALUE TEXT NOT NULL,
            TYPE TEXT NOT NUll,
            FOREIGN KEY(MATKEY) REFERENCES MATERIAL(PKEY)
        );
        """,
        "CREATE INDEX IF NOT EXISTS MFCZ_INDEX ON MASSFRACTION(MFZ, MFC);",
    )
    for cmd in commands
        SQLite.execute(db, cmd)
    end
end

"""
    Base.write(db::SQLite.DB, ::Type{Material}, mat::Material)::Int

Add a `Material` to the SQLite database.  Will not overwrite a previously define definition.
To replace a definition, first `delete!(db, Material, matname)`.
Returns the database key associated with the `Material`.
"""
function Base.write(db::SQLite.DB, ::Type{Material}, mat::Material)::Int
    haskey(db, Material, mat.name) && error("$(mat.name) has already been defined as $(mat).")
    SQLite.transaction(db) do
        stmt1 = SQLite.Stmt(db, "INSERT INTO MATERIAL (NAME, DENSITY) VALUES ( ?, ? );")
        r = DBInterface.execute(stmt1, (name(mat), get(mat, :Density, missing)))
        pkey = DBInterface.lastrowid(r)
        stmt3 = SQLite.Stmt(db, "INSERT INTO MASSFRACTION ( MATKEY, MFZ, MFC, MFUC, MFA ) VALUES ( ?, ?, ?, ?, ? );")
        foreach(elm->DBInterface.execute(stmt3, (pkey, z(elm), value(mat[elm]), σ(mat[elm]), get(mat.a, elm, 0.0))), keys(mat))
        stmt4 = SQLite.Stmt(db, "INSERT INTO MATERIALPROPERTY ( MATKEY, KEY, VALUE, TYPE ) VALUES ( ?, ?, ?, ? );")
        for (key, val) in filter(p->p[2] isa AbstractString, mat.properties)
            DBInterface.execute(stmt4, (pkey, repr(key)[2:end], val, "String"))
        end
        return pkey
    end
end

"""
    Base.read(db::SQLite.DB, ::Type{Material}, pkey::Int)::Material
    Base.read(db::SQLite.DB, ::Type{Material}, matname::AbstractString)::Material

Read a `Material` from the database by index or by name.
"""
function Base.read(db::SQLite.DB, ::Type{Material}, pkey::Int)::Material
    stmt1 = SQLite.Stmt(db, "SELECT * FROM MATERIAL WHERE PKEY=?;")
    q1 = DBInterface.execute(stmt1, (pkey, ))
    SQLite.done(q1) && error("No known material with pkey = '$(pkey)'.")
    r1 = first(q1)
    row, den = r1[:PKEY], r1[:DENSITY]
    stmt2 = SQLite.Stmt(db, "SELECT * FROM MASSFRACTION WHERE MATKEY=?;")
    massfrac, aa = Dict{Element,UncertainValue}(), Dict{Element,Float64}()
    for r2 in DBInterface.execute(stmt2, (row, ))
        @assert r2[:MATKEY]==row
        z, c, uc, a = elements[r2[:MFZ]], r2[:MFC], r2[:MFUC], r2[:MFA]
        massfrac[z] = uv(c,uc)
        a>0.0 && (aa[z]=a)
    end
    props=Dict{Symbol,Any}()
    (!ismissing(den)) && (props[:Density]=den)
    stmt3 = SQLite.Stmt(db, "SELECT * FROM MATERIALPROPERTY WHERE MATKEY=?;")
    for r3 in DBInterface.execute(stmt3, (row, ))
        @assert r3[:MATKEY]==row
        key, val = r3[:KEY], r3[:VALUE]
        props[Symbol(key)]=val
    end
    return Material(r1[:NAME], massfrac, aa, props)
end
function Base.read(db::SQLite.DB, ::Type{Material}, matname::AbstractString)::Material
    stmt1 = SQLite.Stmt(db, "SELECT * FROM MATERIAL WHERE NAME=?;")
    q1 = DBInterface.execute(stmt1, (matname, ))
    SQLite.done(q1) && error("There is no material named $matname in the database.")
    r1 = first(q1)
    row, den = r1[:PKEY], r1[:DENSITY]
    stmt2 = SQLite.Stmt(db, "SELECT * FROM MASSFRACTION WHERE MATKEY=?;")
    massfrac, aa = Dict{Element,UncertainValue}(), Dict{Element,Float64}()
    for r2 in DBInterface.execute(stmt2, (row, ))
        @assert r2[:MATKEY]==row
        z, c, uc, a = elements[r2[:MFZ]], r2[:MFC], r2[:MFUC], r2[:MFA]
        massfrac[z] = uv(c,uc)
        a>0.0 && (aa[z]=a)
    end
    props=Dict{Symbol,Any}()
    (!ismissing(den)) && (props[:Density]=den)
    stmt3 = SQLite.Stmt(db, "SELECT * FROM MATERIALPROPERTY WHERE MATKEY=?;")
    for r3 in DBInterface.execute(stmt3, (row, ))
        @assert r3[:MATKEY]==row
        key, val = r3[:KEY], r3[:VALUE]
        props[Symbol(key)]=val
    end
    return Material(r1[:NAME], massfrac, aa, props)
end

"""
    Base.haskey(db::SQLite.DB, ::Type{Material}, matname::AbstractString)::Bool

Is a `Material` named `matname` in the database?
"""
function Base.haskey(db::SQLite.DB, ::Type{Material}, matname::AbstractString)::Bool
    stmt1 = SQLite.Stmt(db, "SELECT PKEY FROM MATERIAL WHERE NAME=?;")
    q1 = DBInterface.execute(stmt1, (matname, ))
    return !SQLite.done(q1)
end

NeXLCore.material(db::SQLite.DB, matname::AbstractString) =
    read(db, Material, matname)

"""
    Base.delete!(db::SQLite.DB, ::Type{Material}, matname::AbstractString)

Delete the named `Material` from the database.
"""
function Base.delete!(db::SQLite.DB, ::Type{Material}, matname::AbstractString)
    DBInterface.transaction(db) do
        stmt1 = SQLite.Stmt(db, "SELECT PKEY FROM MATERIAL WHERE NAME=?;")
        q1 = DBInterface.execute(stmt1, (matname, ))
        for r1 in q1
            SQLite.transaction(db) do
                stmt1 = SQLite.Stmt(db, "DELETE FROM MASSFRACTION where MATKEY=?;")
                DBInterface.execute(stmt1, (r1[:PKEY], ))
                stmt2 = SQLite.Stmt(db, "DELETE FROM MATERIALPROPERTY where MATKEY=?;")
                DBInterface.execute(stmt2, (r1[:PKEY], ))
                stmt3 = SQLite.Stmt(db, "DELETE FROM MATERIAL where PKEY=?;")
                DBInterface.execute(stmt3, (r1[:PKEY], ))
            end
        end
    end
end

"""
    Base.findall(db::SQLite.DB, ::Type{Material}, prs::Pair{Element, <:Tuple{<:Real,<:Real}}...)::Vector{<:Material}
    Base.findall(db::SQLite.DB, ::Type{Material}, filt::Dict{Element, <:Tuple{<:Real,<:Real}})::Vector{<:Material}
    Base.findall(db::SQLite.DB, mat::Material, tol::Float64)::Vector{<:Material}

Returns all the `Material`s in the database that have mass fractions for the specified elements within a range. The
first two versions take a collection of `Element`, `Tuple` pairs that specify the range of mass fractions for each
element.

Examples:

    julia> findall(db, Material, Dict(n"Bi"=>(0.1,0.5), n"Cu"=>(0.2,0.3)))
    3-element Vector{Material{UncertainValue, Float64}}:
     Favreauite[Bi=0.1488,Se=0.2249,Cu=0.2715,Pb=0.1475,O=0.2051,H=0.0022]
     Miharaite[Bi=0.2275,Cu=0.2767,S=0.2094,Fe=0.0608,Pb=0.2256]
     Mrázekite[Bi=0.4641,Cu=0.2117,P=0.0688,O=0.2487,H=0.0067]

    julia> findall(db, mat"NaAlSi3O8", 0.01)
    4-element Vector{Material{UncertainValue, Float64}}:
     Albite[Al=0.1029,Na=0.0877,Si=0.3213,O=0.4881]
     Kumdykolite[Al=0.1029,Na=0.0877,Si=0.3213,O=0.4881]
     Lingunite[Al=0.1029,Na=0.0877,Si=0.3213,O=0.4881]
     Monalbite[Al=0.1029,Na=0.0877,Si=0.3213,O=0.4881]
"""
Base.findall(db::SQLite.DB, ::Type{Material}, prs::Pair{Element, <:Tuple{<:Real,<:Real}}...)::Vector{<:Material} =
    findall(db, Material, Dict(prs...))

function Base.findall(db::SQLite.DB, ::Type{Material}, filt::Dict{Element, <:Tuple{<:Real,<:Real}})::Vector{<:Material}
    cmds, args = String[], Any[]
    for (elm, ci) in filt
        push!(cmds, "SELECT MATKEY FROM MASSFRACTION WHERE MFZ=? AND MFC>=? and MFC<=?")
        append!(args, ( z(elm), ci[1], ci[2] ))
    end
    stmt = SQLite.Stmt(db, join(cmds," INTERSECT ")*";")
    q = DBInterface.execute(stmt, args)
    return [ read(db, Material, r[:MATKEY]) for r in q ]
end
function Base.findall(db::SQLite.DB, mat::Material, tol::Float64)::Vector{<:Material}
    interval(el) = (max(0.0,value(mat[el])-tol), value(mat[el])+tol)
    return findall(db, Material, Dict(elm=>interval(elm) for elm in keys(mat)))
end

"""
    Base.findall(db::SQLite.DB, ::Type{Material}, like::AbstractString)::Vector{String}

Returns the names of all `Material`s in the database which match the `like` string.

`like` uses SQL like syntax where `%` and `_` have special meanings. Searches are case-insenstive.

  * '%' matches zero or more characters
  * '_' matches one character

Using the RUFF database:

  * `findall(db, Material, "Al%")` matches all 37 materials starting with the letters 'a' then 'l' (case-insenstive).
  * `findall(db, Material, "Al%ite")` matches all materials starting with the letters 'a' then 'l' and ending with 'ite'.
  * `findall(db, Material, "Al__rsite")` matches "Alforsite" and "Alpersite" while `findall(db, Material, "Al%rsite")` matches these two plus `Alarsite`.
"""
function Base.findall(db::SQLite.DB, ::Type{Material}, like::AbstractString)::Vector{String}
    s=SQLite.Stmt(db, "SELECT NAME FROM MATERIAL WHERE NAME LIKE ?;")
    map(r->r[1], DBInterface.execute(s, ( like, )))
end

# md = loadmineraldata(true)
# foreach(r->write(db,Material,r[:Material]),filter(r->r[:Material] isa Material, eachrow(md)))


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Libuuid"
version = v"2.34"

# Collection of sources required to build FriBidi
sources = [
    ArchiveSource("https://mirrors.edge.kernel.org/pub/linux/utils/util-linux/v$(version.major).$(version.minor)/util-linux-$(version.major).$(version.minor).tar.xz",
                  "743f9d0c7252b6db246b659c1e1ce0bd45d8d4508b4dfa427bbb4a3e9b9f62b5"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/util-linux-*/
./configure --prefix=$prefix --host=$target --disable-all-programs --enable-libuuid
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [p for p in supported_platforms() if !(p isa Union{Windows,MacOS})]

# The products that we will ensure are always built
products = [
    LibraryProduct("libuuid", :libuuid)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)