text:
@testset "AD: Others" begin
    if GROUP == "All" || GROUP == "Tracker"
        @testset "TuringUniform" begin
            @test logpdf(TuringUniform(), param(0.5)) == 0
        end

        @testset "Semicircle" begin
            @test Tracker.data(logpdf(Semicircle(1.0), param(0.5))) == logpdf(Semicircle(1.0), 0.5)
        end
    end

    @testset "logsumexp" begin
        x = rand(3)
        test_reverse_mode_ad(logsumexp, randn(), x; rtol=1e-8, atol=1e-6)
    end

    @testset "zygote_ldiv" begin
        A = to_posdef(rand(3, 3))
        B = to_posdef(rand(3, 3))

        test_reverse_mode_ad(randn(3, 3), A, B) do A, B
            return DistributionsAD.zygote_ldiv(A, B)
        end
    end

    @testset "logdet" begin
        N = 7
        B = randn(N, N)

        test_reverse_mode_ad(randn(), B; rtol=1e-8, atol=1e-6) do B
            return logdet(cholesky(to_posdef(B)))
        end
        test_reverse_mode_ad(randn(), B; rtol=1e-8, atol=1e-6) do B
            return logdet(cholesky(Symmetric(to_posdef(B))))
        end
    end

    @testset "fill" begin
        if GROUP == "All" || GROUP == "Tracker"
            @test fill(param(1.0), 3) isa TrackedArray
        end

        test_reverse_mode_ad(x->fill(x, 7), randn(7), randn())
        test_reverse_mode_ad(x->fill(x, 7, 11), randn(7, 11), randn())
        test_reverse_mode_ad(x->fill(x, 7, 11, 13), rand(7, 11, 13), randn())
    end

    @testset "Tracker, Zygote and ReverseDiff + MvNormal" begin
        N = 7
        m = rand(N)
        B = randn(N, N)
        x = rand(TuringDenseMvNormal(m, to_posdef(B)))

        test_reverse_mode_ad(randn(), m, B, x) do m, B, x
            return logpdf(MvNormal(m, to_posdef(B)), x)
        end
        test_reverse_mode_ad(randn(), m, B, x) do m, B, x
            return logpdf(TuringMvNormal(m, to_posdef(B)), x)
        end
        test_reverse_mode_ad(randn(), m, B, x) do m, B, x
            return logpdf(TuringMvNormal(m, Symmetric(to_posdef(B))), x)
        end
    end

    @testset "adapt_randn" begin
        rng = MersenneTwister()
        n = 50
        dims = (10, 30)
        for T in (Float32, Float64)
            if GROUP == "All" || GROUP == "ForwardDiff"
                let
                    x = [ForwardDiff.Dual(rand(rng, T)) for _ in 1:n]
                    test_adapt_randn(rng, x, T, dims...)
                end
            end
            if GROUP == "All" || GROUP == "Tracker"
                let
                    x = Tracker.TrackedArray(rand(rng, T, 50))
                    test_adapt_randn(rng, x, T, dims...)
                end
            end
            if GROUP == "All" || GROUP == "ReverseDiff"
                let
                    v = rand(rng, T, n)
                    d = rand(Int, n)
                    tp = ReverseDiff.InstructionTape()
                    x = ReverseDiff.TrackedArray(v, d, tp)
                    test_adapt_randn(rng, x, T, dims...)
                end
            end
        end
    end
end


text:
function minmeanmedianmax(vec)
    (minimum(vec),round(mean(vec),2),median(vec),maximum(vec))
end
function statcols(data)
    num_col = size(data,2)
    for ind in 1:num_col
        tmp = data[:,ind]
        print("min mean median max for $l class: ")
        println(minmeanmedianmax(tmp))
    end
end

ext: jl


text:
module PATH

import MathOptInterface
import SparseArrays

include(joinpath(dirname(@__DIR__), "deps", "deps.jl"))

include("C_API.jl")
include("MOI_wrapper.jl")

end


text:
using Base.Test

using RaggedData:
    # old interface
    RaggedCounter, collate_index_keys,
    RaggedCollate, next_index!,
    RaggedIndex, RaggedColumns, RaggedColumn, _subset,
    # experimental interface
    ordered_counts, contiguous_ranges, contiguous_invperm, contiguous_invperm!

using DataStructures: OrderedDict

@testset "ragged data processing" begin
    # initialize and countr
    rc = RaggedCounter(Int,Int32)
    @test rc(3) == 3
    @test rc(2) == 2
    @test rc(1) == 1
    for _ in 1:3
        @test rc(1) == 1
    end
    for _ in 1:10
        @test rc(2) == 2
    end
    for _ in 1:7
        @test rc(3) == 3
    end
    @test length(rc) == 3
    @test count(rc) == 3+3+10+7
    # test index and collation
    coll, ix, ks = collate_index_keys(rc, true)
    @test ks == Int32[3, 2, 1]
    # index
    @test length(ix) ≡ 3
    @test size(ix) ≡ (3,)
    @test eltype(ix) ≡ UnitRange{Int32}
    @test ix[1] ≡ UnitRange{Int32}(1, 8)
    @test ix[2] ≡ UnitRange{Int32}(9, 19)
    @test ix[3] ≡ UnitRange{Int32}(20, 23)
    @test indices(ix) ≡ indices(ones(3))
    # collation
    @test next_index!(coll, 3) == 1
    for i in 9:15
        @test next_index!(coll, 2) == i
    end
    for i in 20:23
        @test next_index!(coll, 1) == i
    end
    for i in 16:19
        @test next_index!(coll, 2) == i
    end
    for i in 2:8
        @test next_index!(coll, 3) == i
    end
    # equality testing
    @test isequal(rc, deepcopy(rc))
    @test isequal(coll, deepcopy(coll))
    @test isequal(ix, deepcopy(ix))
end

@testset "ragged index subset calculations" begin
    ix = RaggedIndex(cumsum([5,3,2,8]))
    sub_ix, sub_I = _subset(ix, 3:4)
    @test sub_ix == RaggedIndex(cumsum([2,8]))
    @test sub_I == collect(9:18)
end

@testset "ragged columns" begin
    ix = RaggedIndex(cumsum([5,3,2,7]))
    v1 = collect(1:17)
    v2 = collect(StepRangeLen(18.0, 1.0, 17))

    # columns
    rc = RaggedColumns(ix, (v1, v2))
    @test length(rc) == 4
    @test count(rc) == 17
    @test eltype(rc) == typeof(((@view v1[1:1]), (@view v2[1:1]))) == typeof(rc[1])
    @test size(rc) == (4,)
    @test size(rc, 1) == 4
    for i in indices(rc, 1)
        j = ix[i]
        @test rc[i] == (v1[j], v2[j])
    end
    @test rc[3:4] == RaggedColumns(RaggedIndex(cumsum([2,7])), (v1[9:17], v2[9:17]))
    @test collect(rc) ==
        [([1, 2, 3, 4, 5], [18.0, 19.0, 20.0, 21.0, 22.0]),
         ([6, 7, 8], [23.0, 24.0, 25.0]),
         ([9, 10], [26.0, 27.0]),
         ([11, 12, 13, 14, 15, 16, 17], [28.0, 29.0, 30.0, 31.0, 32.0, 33.0, 34.0])]

    # subset of columns
    rc1 = RaggedColumns(rc, [1])
    @test length(rc1) == 4
    @test count(rc1) == 17
    @test eltype(rc1) == typeof(((@view v1[1:1]),)) == typeof(rc1[1])
    @test size(rc1) == (4,)
    @test size(rc1, 1) == 4
    for i in indices(rc1, 1)
        j = ix[i]
        @test rc1[i] == (v1[j],)
    end
    @test_throws ArgumentError RaggedColumns(rc, 1)
    @test collect(rc1) ==
        [([1, 2, 3, 4, 5], ),
         ([6, 7, 8], ),
         ([9, 10], ),
         ([11, 12, 13, 14, 15, 16, 17], )]

    # single column
    rc2 = RaggedColumn(rc, 1)
    @test length(rc2) == 4
    @test count(rc2) == 17
    @test eltype(rc2) == typeof(@view v1[1:1]) == typeof(rc2[1])
    @test size(rc2) == (4,)
    @test size(rc2, 1) == 4
    for i in indices(rc2, 1)
        j = ix[i]
        @test rc2[i] == v1[j]
    end
    @test_throws MethodError RaggedColumn(rc, 1:1)
    @test collect(rc2) ==
        Any[[1, 2, 3, 4, 5], [6, 7, 8], [9, 10], [11, 12, 13, 14, 15, 16, 17]]
    @test rc2[2:3] == [[6, 7, 8], [9, 10]]
end

@testset "contiguous permutations" begin
    for _ in 1:30
        N = rand(5:20)          # number of distinct elements
        K = rand(100:200)       # number of additional observations
        T = rand([Int64, Int32, Int16, UInt64, UInt32, UInt16])
        xs = rand(1:N, N + K)
        xs[1:N] .= 1:N              # to establish order
        c = OrderedDict{Int, T}()
        for x in xs
            c[x] = get(c, x, 0) + 1
        end
        c2 = ordered_counts(xs)
        @test c2 == c
        rs = contiguous_ranges(c)
        @test eltype(rs) ≡ UnitRange{T}
        for (i, r) in enumerate(rs)
            @test length(r) == sum(xs .== i)
            if i > 1
                @test rs[i-1].stop + 1 == r.start
            else
                @test r.start == 1
            end
        end
        p = contiguous_invperm(xs, c)
        @test length(p) == length(xs)
        @test eltype(p) ≡ T
        @test sort(p) == collect(indices(xs, 1))
        z = collect(enumerate(xs))
        zp = similar(z)
        zp[p] .= z
        for (i, r) in enumerate(rs)
            @test issorted(first.(zp[r]))
            @test all(last.(zp[r]) .== i)
        end
    end
end


text:
module DepositionMTK
import Base.show, StaticArrays
export GasData, r_s, r_dc, r_mx, r_smx, r_lux, r_clx, r_gsx, SurfaceResistance, inf, r_i, r_lu, r_ac, r_gsS, r_gsO, r_clS, r_clO, So2Data, O3Data, No2Data, NoData, Hno3Data, H2o2Data, GasData, AldData, HchoData, OpData, PaaData, OraData, Nh3Data, PanData, Hno2Data

	using StaticArrays

	const inf = 1.e25

	# r_i represents the minimum bulk canopy resistances for water vapor.
	const r_i = SA_F32[
	inf 60 120 70 130 100 inf inf 80 100 150
	inf inf inf inf 250 500 inf inf inf inf inf
	inf inf inf inf 250 500 inf inf inf inf inf
	inf inf inf inf 400 800 inf inf inf inf inf
	inf 120 240 140 250 190 inf inf 160 200 300]

	# r_lu signifies leaf cuticles in healthy vegetation and otherwise the outer surfaces in the upper canopy.
	const r_lu = SA_F32[
	inf 2000 2000 2000 2000 2000 inf inf 2500 2000 4000
	inf 9000 9000 9000 4000 8000 inf inf 9000 9000 9000
	inf inf 9000 9000 4000 8000 inf inf 9000 9000 9000
	inf inf inf inf 6000 9000 inf inf 9000 9000 9000
	inf 4000 4000 4000 2000 3000 inf inf 4000 4000 8000]

	# r_ac signifies transfer that depends only on canopy height and density.
	const r_ac = SA_F32[
	100.0 200 100 2000 2000 2000 0 0 300 150 200
	100 150 100 1500 2000 1700 0 0 200 120 140
	100 10 100 1000 2000 1500 0 0 100 50 120
	100 10 10 1000 2000 1500 0 0 50 10 50
	100 50 80 1200 2000 1500 0 0 200 60 120]

	# r_gs signifies uptake at the "ground" by soil, leaf litter, snow, water etc. 'S' and 'O' stand for SO2 and O3 respectively.
	const r_gsS = SA_F32[
	400.0 150 350 500 500 100 0 1000 0 220 40
	400 200 350 500 500 100 0 1000 0 300 400
	400 150 350 500 500 200 0 1000 0 200 400
	100 100 100 100 100 100 0 1000 100 100 50
	500 150 350 500 500 200 0 1000 0 250 40]

	const r_gsO = SA_F32[
	300.0 150 200 200 200 300 2000 400 1000 180 200
	300 150 200 200 200 300 2000 400 800 180 200
	300 150 200 200 200 300 2000 400 1000 180 20
	600 3500 3500 3500 3500 3500 2000 400 3500 3500 3500
	300 150 200 200 200 300 2000 400 1000 180 200]

	# r_cl is meant to account for uptake pathways at the leaves, bark, etc. 'S' and 'O' stand for SO2 and O3 respectively.
	const r_clS = SA_F32[
	inf 2000 2000 2000 2000 2000 inf inf 2500 2000 4000
	inf 9000 9000 9000 2000 4000 inf inf 9000 9000 9000
	inf inf 9000 9000 3000 6000 inf inf 9000 9000 9000
	inf inf inf 9000 200 400 inf inf 9000 inf 9000
	inf 4000 4000 4000 2000 3000 inf inf 4000 4000 8000]

	const r_clO = SA_F32[
	inf 1000 1000 1000 1000 1000 inf inf 1000 1000 1000
	inf 400 400 400 1000 600 inf inf 400 400 400
	inf 1000 400 400 1000 600 inf inf 800 600 600
	inf 1000 1000 400 1500 600 inf inf 800 1000 800
	inf 1000 500 500 1500 700 inf inf 600 800 800]

	# Holder for gas properties from Wesely (1989) Table 2.'
	struct GasData
		Dh2oPerDx::AbstractFloat
		Hstar::AbstractFloat
		Fo::AbstractFloat
	end

	const So2Data = GasData(1.9, 1.e5, 0)
	const O3Data  = GasData(1.6, 0.01, 1)
	const No2Data = GasData(1.6, 0.01, 0.1) # Wesely (1989) suggests that,
	# in general, the sum of NO and NO2 should be considered rather
	# than NO2 alone because rapid in-air chemical reactions can cause
	# a significant change of NO and NO2 vertical fluxes between the
	# surface and the point at which deposition velocities are applied,
	# but the sum of NO and NO2 fluxes should be practically unchanged.
	const NoData   = GasData(1.3, 3.e-3, 0) # Changed according to Walmsley (1996)
	const Hno3Data = GasData(1.9, 1.e14, 0)
	const H2o2Data = GasData(1.4, 1.e5, 1)
	const AldData  = GasData(1.6, 15, 0)     # Acetaldehyde (aldehyde class)
	const HchoData = GasData(1.3, 6.e3, 0)   # Formaldehyde
	const OpData   = GasData(1.6, 240, 0.1)
	# Methyl hydroperoxide (organic peroxide class)
	const PaaData  = GasData(2.0, 540, 0.1)  # Peroxyacetyl nitrate
	const OraData  = GasData(1.6, 4.e6, 0)   # Formic acid (organic acid class)
	const Nh3Data  = GasData(0.97, 2.e4, 0)  # Changed according to Walmsley (1996)
	const PanData  = GasData(2.6, 3.6, 0.1)  # Peroxyacetyl nitrate
	const Hno2Data = GasData(1.6, 1.e5, 0.1) # Nitrous acid

	# Calculate bulk canopy stomatal resistance [s m-1] based on Wesely (1989) equation 3 when given the solar irradiation (G [W m-2]), the surface air temperature (Ts [°C]), the season index (iSeason), the land use index (iLandUse), and whether there is currently rain or dew.
	function r_s(G::T, Ts::T, iSeason::Int, iLandUse::Int, rainOrDew::Bool) where T<:AbstractFloat
		rs = 0.0
		if Ts >= 39.9 || Ts <= 0.1
			rs = inf
		else
			rs = r_i[iSeason, iLandUse] * (1 + (200.0 * 1.0 / (G + 1.0))^2) *
				(400.0 * 1.0 / (Ts * (40.0 - Ts)))
		end
		# Adjust for dew and rain (from "Effects of dew and rain" section).
		if rainOrDew
			rs *= 3
		end
		return rs
	end

	# Calculate the resistance from the effects of mixing forced by buoyant convection when sunlight heats the ground or lower canopy and by penetration of wind into canopies on the sides of hills [s m-1] when given the solar irradiation (G [W m-2]) and the slope of the local terrain (θ [radians]). From Wesely (1989) equation 5.
	function r_dc(G::T, θ::T) where T<:AbstractFloat
		return 100.0 * (1.0 + 1000.0/(G+10.0)) / (1.0 + 1000.0*θ)
	end

	#Calculate mesophyll resistance [s m-1] based on Wesely (1989) equation 6 when given the effective Henry's law coefficient (Hstar [M atm-1]) and the reactivity factor (fo [-]), both available in Wesely (1989) table 2.
	function r_mx(Hstar::T, fo::T) where T<:AbstractFloat
		return 1.0 / (Hstar/3000.0 + 100.0*fo)
	end

	#Calculate combined minimum stomatal and mesophyll resistance [s m-1] based on Wesely (1989) equation 4 when given stomatal resistance (r_s [s m-1]), ratio of water to chemical-of-interest diffusivities (Dh2oPerDx [-]), and mesophyll resistance (r_mx [s m-1]).
	function r_smx(r_s::T, Dh2oPerDx::T, r_mx::T) where T<:AbstractFloat
		return r_s * Dh2oPerDx + r_mx
	end

	# Calculate the resistance of the outer surfaces in the upper canopy (leaf cuticular resistance in healthy vegetation) based on Wesely (1989) equations 7 and 10-14 when given the effective Henry's law coefficient (Hstar [M atm-1]), the reactivity factor (fo [-]) (both available in Wesely (1989) table 2), the season index (iSeason), the land use index (iLandUse), whether there is currently rain or dew, and whether the chemical of interest is either SO2 or O3.
	function r_lux(Hstar::T, fo::T, iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) where T<:AbstractFloat
		rlux = 0.0
		if dew && (iSeason != 4) # Dew doesn't have any effect in the winter
			if isSO2
				if iLandUse == 1
					rlux = 50.0 # equation 13 and a half
				else
					rlux = 100 # equation 10.
				end
			elseif isO3
				# equation 11
				rlux = 1.0 / (1.0/3000.0 + 1.0/(3*r_lu[iSeason, iLandUse]))
			else
				rluO = 1.0 / (1.0/3000.0 + 1.0/(3*r_lu[iSeason, iLandUse])) #equation 11
				rlux = 1.0 / (1.0/(3*r_lu[iSeason, iLandUse]/(1.0e-5*Hstar+fo))
					+ 1.0e-7*Hstar + fo/rluO) # equation 14, modified to match Walmsley eq. 5g
			end
		elseif rain && (iSeason != 4)
			if isSO2
				if iLandUse == 1
					rlux = 50 #equation 13 and a half
				else
					# equation 12
					rlux = 1.0 / (1.0/5000.0 + 1.0/(3*r_lu[iSeason, iLandUse]))
				end
			elseif isO3
				# equation 13
				rlux = 1.0 / (1.0/1000.0 + 1.0/(3*r_lu[iSeason, iLandUse]))
			else
				rluO = 1.0 / (1.0/1000.0 + 1.0/(3*r_lu[iSeason, iLandUse])) # equation 13
				rlux = 1.0 / (1.0/(3*r_lu[iSeason, iLandUse]/(1.e-5*Hstar+fo)) + 1.0e-7*Hstar +fo/rluO) # equation 14, modified to match Walmsley eq. 5g
			end
		else
			rlux = r_lu[iSeason, iLandUse] / (1.0e-5*Hstar + fo)
		end
		return rlux
	end

	# Calculate the resistance of the exposed surfaces in the lower portions of structures (canopies, buildings) above the ground based on Wesely (1989) equation 8 when given the effective Henry's law coefficient (Hstar [M atm-1]), the reactivity factor (fo [-]) (both available in Wesely (1989) table 2), the season index (iSeason), and the land use index (iLandUse).
	function r_clx(Hstar::T, fo::T, iSeason::Int, iLandUse::Int) where T<:AbstractFloat
		return 1.0 / (Hstar/(1.0e5*r_clS[iSeason, iLandUse]) +
			fo/r_clO[iSeason, iLandUse])
	end

	# Calculate the resistance to uptake at the 'ground' surface based on Wesely (1989) equation 9 when given the effective Henry's law coefficient (Hstar [M atm-1]), the reactivity factor (fo [-]) (both available in Wesely (1989) table 2), the season index (iSeason), and the land use index (iLandUse).
	function r_gsx(Hstar::T, fo::T, iSeason::Int, iLandUse::Int) where T<:AbstractFloat
		return 1.0 / (Hstar/(1.0e5*r_gsS[iSeason, iLandUse]) +
			fo/r_gsO[iSeason, iLandUse])
	end

	# Calculates surface resistance to dry depostion [s m-1] based on Wesely (1989) equation 2 when given information on the chemical of interest (gasData), solar irradiation (G [W m-2]), the surface air temperature (Ts [°C]), the slope of the local terrain (Θ [radians]), the season index (iSeason), the land use index (iLandUse), whether there is currently rain or dew, and whether the chemical of interest is either SO2 (isSO2) or O3 (isO3).
	# From Wesely (1989) regarding rain and dew inputs: "A direct computation of the surface wetness would be most desirable, e.g. by estimating the amount of free surface water accumulated and then evaporated. Alternatively, surface relative humidity might be a useful	index. After dewfall and rainfall events are completed, surface wetness	often disappears as a result of evaporation after approximately 2	hours of good atmospheric mixing, the period of time recommended earlier (Sheih et al., 1986)".
	function SurfaceResistance(gasData::GasData, G::T, Ts::T, θ::T,
		iSeason::Int, iLandUse::Int, rain::Bool, dew::Bool, isSO2::Bool, isO3::Bool) where T<:AbstractFloat
		rs = r_s(G, Ts, iSeason, iLandUse, rain || dew)
		rmx = r_mx(gasData.Hstar, gasData.Fo)
		rsmx = r_smx(rs, gasData.Dh2oPerDx, rmx)
		rdc = r_dc(G, θ)
		rlux = r_lux(gasData.Hstar, gasData.Fo, iSeason, iLandUse,
			rain, dew, isSO2, isO3)
		rclx = 0.0
		rgsx = 0.0
		if isSO2
			rclx = r_clS[iSeason, iLandUse]
			rgsx = r_gsS[iSeason, iLandUse]
		elseif isO3
			rclx = r_clO[iSeason, iLandUse]
			rgsx = r_gsO[iSeason, iLandUse]
		else
			rclx = r_clx(gasData.Hstar, gasData.Fo, iSeason, iLandUse)
			rgsx = r_gsx(gasData.Hstar, gasData.Fo, iSeason, iLandUse)
		end

		rac = r_ac[iSeason, iLandUse]

		# Correction for cold temperatures from page 4 column 1.
		if Ts < 0.0
			correction = 1000.0 * exp(-Ts-4) # [s m-1] #mark
			rlux += correction
			rclx += correction
			rgsx += correction
		end

		r_c = 1.0 / (1.0/(rsmx) + 1.0/rlux + 1.0/(rdc+rclx) + 1.0/(rac+rgsx))
		r_c = max(r_c, 10.0) # From "Results and conclusions" section to avoid extremely high deposition velocities over extremely rough surfaces.
		r_c = min(r_c, 9999.0)
		return r_c
	end
end


text:
# S is X,Y
# T is maxpasses
# p: # of features
# n: # of data points
# k: size of minibatch
function pegasos(X::SparseOrFullMat,
	             Y::Vector;
	             k::Integer = 5,
	             lambda::Real = 0.1,
	             T::Integer = 100)
	# p features, n observations
	p, n = size(X)

	# Initialize weights so norm(w) <= 1 / sqrt(lambda)
	w = randn(p)
	sqrtlambda = sqrt(lambda)
	normalizer = sqrtlambda * norm(w)
	for j in 1:p
		w[j] /= normalizer
	end

	# Allocate storage for repeated used arrays
	deltaw = Array(Float64, p)
	w_tmp = Array(Float64, p)

	# Loop
	for t in 1:T
		# Calculate stepsize parameters
		alpha = 1.0 / t
		eta_t = 1.0 / (lambda * t)

		# Calculate scaled sum over misclassified examples
		# Subgradient over minibatch of size k
		fill!(deltaw, 0.0)
		for i in 1:k
			# Select a random item from X
			# This is one element of At of S
			index = rand(1:n)

			# Test if prediction isn't sufficiently good
			# If so, current item is element of At+
			pred = Y[index] * inner_prod(w, X[:, index])
			if pred < 1.0
				# Update subgradient
				add_vec!(deltaw, X[:, index], e -> Y[index] * e)
			end
		end

		# Rescale subgradient
		for j in 1:p
			deltaw[j] *= (eta_t / k)
		end

		# Calculate tentative weight-update
		for j in 1:p
			w_tmp[j] = (1.0 - alpha) * w[j] + deltaw[j]
		end

		# Find projection of weights into L2 ball
		proj = min(1.0, 1.0 / (sqrtlambda * norm(w_tmp)))
		for j in 1:p
			w[j] = proj * w_tmp[j]
		end
	end

	return SVMFit(w, T, true)
end


text:
export laplacian_loss, edge_loss

"""
    laplacian_loss(m::TriMesh)

Computes the laplacian smoothning metrics for TriMesh `m`

"""
function laplacian_loss(m::TriMesh)
    L = @ignore get_laplacian_packed(m)
    verts = get_verts_packed(m)
    L = L * cpu(transpose(verts))
    L = _norm(L; dims = 2)
    return mean(L)
end

"""
    edge_loss(m::TriMesh, target_length::Number = 0.0)

Computes mean edge length in TriMesh.

`target_length` is the optional arguments for computing edge length.
"""
function edge_loss(m::TriMesh, target_length::Number = 0.0)
    verts = get_verts_packed(m)
    edges = @ignore get_edges_packed(m)
    v1 = verts[:, edges[:, 1]]
    v2 = verts[:, edges[:, 2]]
    el = (_norm(v1 - v2; dims = 1) .- Float32(target_length)) .^ 2
    loss = mean(el)
    return loss
end

function chamfer_distance(
    m1::TriMesh,
    m2::TriMesh,
    num_samples::Int = 5000;
    w1::Number = 1.0,
    w2::Number = 1.0,
)
    A = sample_points(m1, num_samples)
    B = sample_points(m2, num_samples)
    return _chamfer_distance(A, B, Float32(w1), Float32(w2))
end


text:
module SDM

using LinearAlgebra
using Printf

include("utility.jl")
include("graphics.jl")

export sdm

function sdm(fun, fun_derivative, starting_point, step_size_params, epsilon; display_g=true)
    type = get(step_size_params, "type", "constant")
    method = get(SDM_step_size_methods, type, constant_step_size_sdm)
    return method(fun, fun_derivative, starting_point, step_size_params, epsilon, display_g)
end

function constant_step_size_sdm(fun, fun_derivative, starting_point, step_size_params, epsilon, display_g)
    @printf("Press enter to continue with iterations\n")
    @printf("step\t\tposition\t\tgradient_norm\n")

    gradient = fun_derivative(starting_point[1], starting_point[2])
    x = [starting_point[1], starting_point[2]]

    step_counter = 1
    while(true)

        #we move in the opposite direction of the gradient
        d = -gradient
        #step size is constant
        alpha = get(step_size_params, "step_size", 0.1)
        #move to the next point
        movement = alpha*d
        if display_g
            display(add_vector(x, movement))
        end
        x = x + movement

        gradient_norm = norm(gradient, 2)
        @printf("%d\t\t[%g,%g]\t\t%g", step_counter, x[1], x[2], gradient_norm)

        #we stop whenever the norm of the gradient is lower than a certain epsilon
        if (gradient_norm < epsilon)
            break
        end

        wait_for_key("")

        step_counter += 1

        #calculate the gradient for the new iteration
        gradient = fun_derivative(x[1], x[2])
    end

    @printf("\n")
    return x
end

function quadratic_step_size_sdm(fun, fun_derivative, starting_point, step_size_params, epsilon, display_g)
    @printf("Press enter to continue with iterations\n")
    @printf("step\t\tposition\t\tgradient_norm\n")

    Q = get(step_size_params, "Q", false)
    if Q==false error("Q not set in step_size_params"); return; end;

    gradient = fun_derivative(starting_point[1], starting_point[2])
    x = [starting_point[1], starting_point[2]]

    step_counter = 1
    while(true)

        #we move in the opposite direction of the gradient
        d = -gradient
        #best step size, if the function is a quadratic function
        alpha = (norm(d,2)^2)/(d'*Q*d)
        #move to the next point
        movement = alpha*d
        if display_g
            display(add_vector(x, movement))
        end
        x = x + movement

        gradient_norm = norm(gradient, 2)
        @printf("%d\t\t[%f,%f]\t\t%g", step_counter, x[1], x[2], gradient_norm)

        #we stop whenever the norm of the gradient is lower than a certain epsilon
        if (gradient_norm < epsilon)
            break
        end

        wait_for_key("")

        step_counter += 1

        #calculate the gradient for the new iteration
        gradient = fun_derivative(x[1], x[2])
    end

    @printf("\n")
    return x
end

SDM_step_size_methods = Dict([
    ("constant", constant_step_size_sdm),
    ("quadratic", quadratic_step_size_sdm)
])

end