text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -36.7700820821577 <= q <= 37.6713142442229)

                     add_NL_constraint(m, :(gelu(0.9266287961121966 + 0.9630341445616355*gelu(0.9080159995418975 + 0.10769169268185319*gelu(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.2597058016512084*gelu(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.8597063403720888*gelu(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.17950756837409587*gelu(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + 0.8197021439185934*gelu(-0.9264374759333807 + -0.7028105280156569*gelu(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.6656745963467641*gelu(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.9099166873549671*gelu(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.9082646170966706*gelu(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + -0.8564198456770362*gelu(-0.7118497063916944 + 0.36586202274128254*gelu(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + -0.9129427463398758*gelu(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.5191517354247042*gelu(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.42477634545139953*gelu(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + 0.7868467121279954*gelu(0.8429853377330745 + -0.6377098680254001*gelu(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.04468698484435718*gelu(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.3607403638969684*gelu(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.4004266402265344*gelu(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5])))) + gelu(-0.9209290943615045 + -0.04560727661330999*gelu(0.9080159995418975 + 0.10769169268185319*gelu(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.2597058016512084*gelu(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.8597063403720888*gelu(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.17950756837409587*gelu(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + -0.04187997574232316*gelu(-0.9264374759333807 + -0.7028105280156569*gelu(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.6656745963467641*gelu(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.9099166873549671*gelu(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.9082646170966706*gelu(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + 0.009897273687461627*gelu(-0.7118497063916944 + 0.36586202274128254*gelu(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + -0.9129427463398758*gelu(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.5191517354247042*gelu(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.42477634545139953*gelu(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + -0.4138607422831795*gelu(0.8429853377330745 + -0.6377098680254001*gelu(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.04468698484435718*gelu(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.3607403638969684*gelu(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.4004266402265344*gelu(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5])))) + gelu(0.41055131347137674 + 0.5074173875637427*gelu(0.9080159995418975 + 0.10769169268185319*gelu(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.2597058016512084*gelu(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.8597063403720888*gelu(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.17950756837409587*gelu(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + 0.7865169798440914*gelu(-0.9264374759333807 + -0.7028105280156569*gelu(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.6656745963467641*gelu(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.9099166873549671*gelu(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.9082646170966706*gelu(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + -0.6322955384333278*gelu(-0.7118497063916944 + 0.36586202274128254*gelu(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + -0.9129427463398758*gelu(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.5191517354247042*gelu(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.42477634545139953*gelu(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + 0.03598168048021799*gelu(0.8429853377330745 + -0.6377098680254001*gelu(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.04468698484435718*gelu(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.3607403638969684*gelu(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.4004266402265344*gelu(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5])))) + gelu(-0.7045693885767101 + 0.765322888364711*gelu(0.9080159995418975 + 0.10769169268185319*gelu(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.2597058016512084*gelu(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.8597063403720888*gelu(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.17950756837409587*gelu(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + 0.9087169409161477*gelu(-0.9264374759333807 + -0.7028105280156569*gelu(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.6656745963467641*gelu(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.9099166873549671*gelu(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.9082646170966706*gelu(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + -0.29202574365958966*gelu(-0.7118497063916944 + 0.36586202274128254*gelu(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + -0.9129427463398758*gelu(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.5191517354247042*gelu(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.42477634545139953*gelu(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + -0.14636407330293322*gelu(0.8429853377330745 + -0.6377098680254001*gelu(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.04468698484435718*gelu(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.3607403638969684*gelu(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.4004266402265344*gelu(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Niancat.Games.NiancatGames
using Niancat.Games.NiancatGames: SetPuzzle, Guess, Incorrect, Correct, Rejected, NewPuzzle, GetPuzzle, PuzzleIs, NoPuzzleSet
using Niancat.Games.NiancatGames: CorrectNotification
using Niancat.Languages

@testset "Niancat" begin
    include("guessing_test.jl")
    include("language_test.jl")
    include("setpuzzle_test.jl")
    include("score_test.jl")
    include("state_test.jl")
end

ext: jl


text:
# TODO Clean up the scripts using new commands and function names
# GLOBAL variables in the CanopyRT module have been removed
# Re-organize the tutorial is required
# Tested and passed

# # Canopy Radiative Transfer Example #1
# This example will go through the calculations of hemispheric and bidirectional reflectance as well as emission of Solar Induced Chlorophyll Fluorescence (SIF). The sun sensor geometry is defined by the solar zenith angle SZA, the viewing zenith angle VZA (observer, 0degrees=Nadir looking) and the relative azimuth angle between the sun and observer (0 or 180degrees denoting the principal plane).
#

# First, we include Revise (good for debugging) and Parameters (tools for structures)

#using Revise
using Parameters
#----------------------------------------------------------------------------

# Now include the Land modules

using Land

CanopyRT = Land.CanopyRT

@unpack Canopy4RT,
        CanopyRadiation,
        SoilOpti,
        SolarAngles,
        create_canopy_optical,
        create_leaf_bio,
        create_wl_para_set,
        create_incoming_radiation,
        compute_canopy_geometry!,
        compute_canopy_matrices!,
        computeSIF_Fluxes!,
        simulate_short_wave!,
        derive_canopy_fluxes!,
        compute_thermal_fluxes!,
        fluspect! = CanopyRT
#----------------------------------------------------------------------------

# and unpack some of the mostly used variables or structure (same as writing CanopyRT.leafbio)

#@unpack FT, leafbio, canopy, angles, canOpt, canRad,sunRad,soil, wl, wle, wlf = CanopyRT
#----------------------------------------------------------------------------

# This extracted a number of structures, which are loaded with default values at startup:
# * *leafbio* includes all leaf pigment concentrations, such as Chlorophyll, leaf water, Carotenoids and dry matter content. (mutable struct leafbio)
# * *canopy* includes Canopy properties such as LAI, leaf inclination distribution, clumping factor (mutable struct struct_canopy)
# * *angles* includes SZA, VZA and azimuth (mutable struct struct_angles)
# * *canOpt* saves a variety of canopy optical properties to be computed (mutable struct struct_canopyOptProps)
# * *canRad* includes radiation computed within the Canopy RT module (mutable struct struct_canopyRadiation)
# * *sunRad* provides incoming direct and diffuse incoming solar radiation (mutable struct incomingRadiation)
# * *soil* provides soil albedo and skin temperature
#
# as well as some arrays:
# * *wl* is the wavelength grid (nm) for the solar radiation computation
# * *wle* is the wavelength grid in the photosynthetically active range (a subset of wl, used to define fluorescence excitation and driving photosynthesis)
# * *wlf* is the wavelength grid for SIF (subset of wl). The fluorescence excitation matrices are defined on the basis of wle and wlf.

# Let us just define our own leaf here, we can construct it using standard values but need to provide all dimensions (will be fixed as variables are saved in MArrays):

const FType = Float32

wl_set = create_wl_para_set(FType)
leaf = create_leaf_bio(FType, wl_set.nwl, wl_set.nWlE, wl_set.nWlF);
canopy_rt = Canopy4RT{FType, 20, 3.0}()
canRad_rt = CanopyRadiation{FType, wl_set.nwl, wl_set.nWlF, length(canopy_rt.litab), length(canopy_rt.lazitab), canopy_rt.nlayers}()
canOpt_rt = create_canopy_optical(FType, wl_set.nwl, canopy_rt.nlayers, length(canopy_rt.lazitab), length(canopy_rt.litab); using_marray=false)
sunRad_rt = create_incoming_radiation(FType, wl_set.swl);

## show leaf Chlorophyll content:
@show leaf.Cab
#----------------------------------------------------------------------------

# ---
# ## Leaf Optical Properties
#
# Now we can first run fluspect, which is an extension of the Prospect leaf optical properties program and includes computations of the fluorescence emission responses as well. CanopyRT.optis includes all pigment absorption cross sections as well as refractive index of water to compute leaf optical properties.
# For details, see:
#
# Féret, J.B., Gitelson, A.A., Noble, S.D. and Jacquemoud, S., 2017. PROSPECT-D: Towards modeling leaf optical properties through a complete lifecycle. Remote Sensing of Environment, 193, pp.204-215.
#
# Van der Tol, C., Verhoef, W., Timmermans, J., Verhoef, A. and Su, Z., 2009. An integrated model of soil-canopy spectral radiances, photosynthesis, fluorescence, temperature and energy balance. Biogeosciences, 6(12).
#
# Vilfan, N., Van der Tol, C., Muller, O., Rascher, U. and Verhoef, W., 2016. Fluspect-B: A model for leaf fluorescence, reflectance and transmittance spectra. Remote sensing of environment, 186, pp.596-615.

## Run Fluspect:
CanopyRT.fluspect!(leaf, canopy_rt, wl_set);
#----------------------------------------------------------------------------

# #### Fluorescence excitation matrices
# Now we can visualize the fluorescence excitatiom matrices, which determine how absorbed photons at a specific wavelength *wle* incident on the leaf are emitted backwards (in relation to the incident light but hemispherically diffuse here) or forward. Here, we ignore contributions from Photosystem I and just use a fixed shape as input at the chloroplast level for Photosystem II, hence all spectral shapes are determined by the light distribution within the leaf. Several effects can be observed in the plots below: The first SIF peak at 680nm is much weaker in the forward direction than in the backward direction as more the re-absorption within the leaf is much higher for the forward direction. Similarly, blue incident light penetrates the leaf less than at higher wavelengths, so the chlorophyll re-absorption (backwards) is less re-absorbed by chlorophyll as it emanates from less deep layers.
#

## Define a few wavelengths:
wl_blue = 450.0;
wl_red = 600.0;
wl_FarRed = 740.0;
wl_Red = 685.0;
ind_wle_blue  = argmin(abs.(wl_set.wle .-wl_blue));
ind_wle_red = argmin(abs.(wl_set.wle .-wl_red));
ind_wlf_FR  = argmin(abs.(wl_set.wlf .-wl_FarRed));
ind_wlf_R  = argmin(abs.(wl_set.wlf .-wl_Red));
ind_red = argmin(abs.(wl_set.wl .-wl_Red));
ind_NIR = argmin(abs.(wl_set.wl .-800));

#----------------------------------------------------------------------------

#----------------------------------------------------------------------------

# ---
# ### Leaf reflectance and transmission
# Within Fluspect, leaf reflectance and transmission is calculated (this part is actually identical to what PROSPECT is doing and much less involved than the Fluorescence part). The model basically uses a doubling adding routine to solve for the radiative transfer within the leaf using constituent absorption cross sections and constituent abundances in the leaf (assuming vertically homogenous distribution).
#
# Here, we just create another leaf with different concentrations to show the impact on leaf reflectance and transmission due to changes in chlorophyll content and liquid water content (mostly in the two liquid water absorption bands at 1400 and 1900nm).

## Let's create a leaf with a different Cab and Cw (water) content (Try changing other pigment contents, plot leaf reflectance and transmissions and explain where (spectrally) and why reflectance and transmission changes):
leaf_2 = create_leaf_bio(FType, wl_set.nwl, wl_set.nWlE, wl_set.nWlF);
leaf_2.Cab = 80
leaf_2.Cw = 0.012
## show leaf Chlorophyll content:
## and Run Fluspect:
CanopyRT.fluspect!(leaf_2, canopy_rt, wl_set);

#----------------------------------------------------------------------------

# ---
# ## Moving from the leaf to the entire canopy
#
# Our model here is based on mSCOPE[^1] and we can prescribe different leaf optical properties within a complex canopy. At the moment, the leaf angular distribution is the same vertically (but could be changed in the future).
#
# [^1]: Yang, P., Verhoef, W. and van der Tol, C., 2017. The mSCOPE model: A simple adaptation to the SCOPE model to describe reflectance, fluorescence and photosynthesis of vertically heterogeneous canopies. Remote sensing of environment, 201, pp.1-11.
#

## This is to be changed later but at the moment, we need to generate an Array of leaves, basically for each layer of the canopy

arrayOfLeaves = [create_leaf_bio(FType, wl_set.nwl, wl_set.nWlE, wl_set.nWlF) for i in 1:canopy_rt.nlayers]
for i in 1:canopy_rt.nlayers
    fluspect!(arrayOfLeaves[i], canopy_rt, wl_set)
end

#----------------------------------------------------------------------------

# ### Basic Steps for Canopy RT:
# Below are the basic steps for the canopy radiative transfer (SIF fluxes are dependent on stress levels at each leaf but here we use a standard SIF yield for now):

## Set Soil albedo to 0.2
soil = SoilOpti{FType}(wl_set.wl, FType(0.2)*ones(FType, length(wl_set.wl)), FType[0.1], FType(290.0))
angles = SolarAngles{FType}()

## Compute Canopyoptical properties dependend on sun-sensor and leaf angle distributions:
compute_canopy_geometry!(canopy_rt, angles, canOpt_rt)
## Compute RT matrices with leaf reflectance and transmissions folded in:
compute_canopy_matrices!(arrayOfLeaves, canOpt_rt);
## Perform SW radiation transfer:
simulate_short_wave!(canopy_rt, canOpt_rt, canRad_rt, sunRad_rt, soil);
## Compute outgoing SIF flux (using constant fluorescence efficiency at the chloroplast level)
derive_canopy_fluxes!(canopy_rt, canOpt_rt, canRad_rt, sunRad_rt, soil, arrayOfLeaves, wl_set);
#----------------------------------------------------------------------------

# ### Test a VZA dependence in the principal plane

SIF_FR = Float32[]
SIF_R = Float32[]
reflVIS = Float32[]
reflNIR = Float32[]

## Just running the code over all geometries:

## Set sun SZA to 30 degrees
angles.tts=30
## Set 0 azimuth (principal plane)
angles.psi=0
## LAI of 3:
canopy_rt.LAI = 3
## Define VZA
VZA=collect(-89.5:0.5:89.5)

for VZA_ in VZA
    angles.tto=VZA_
    compute_canopy_geometry!(canopy_rt, angles, canOpt_rt)
    compute_canopy_matrices!(arrayOfLeaves, canOpt_rt);
    simulate_short_wave!(canopy_rt, canOpt_rt, canRad_rt, sunRad_rt, soil);
    computeSIF_Fluxes!(arrayOfLeaves, canOpt_rt, canRad_rt, canopy_rt, soil, wl_set);
    ## Handpicked indices in
    push!(reflVIS, canRad_rt.alb_obs[ind_red])
    push!(reflNIR, canRad_rt.alb_obs[ind_NIR])
    push!(SIF_R , canRad_rt.SIF_obs[ind_wlf_R])
    push!(SIF_FR, canRad_rt.SIF_obs[ind_wlf_FR ])
end
#----------------------------------------------------------------------------

# ---
# ## BRDF sampling
# By going through viewing and azimuth angles, we can construct a full BRDF for reflectance and SIF emissions at different wavelengths:

reflVIS = Float32[]
reflNIR = Float32[]
SIF_FR = Float32[]
SIF_R  = Float32[]
angles.tts=48
angles.psi=0
canopy_rt.LAI=3.22
for psi=0:360
    angles.psi=psi
    for VZA=0:1:85
        angles.tto=VZA

        compute_canopy_geometry!(canopy_rt, angles, canOpt_rt)
        compute_canopy_matrices!(arrayOfLeaves, canOpt_rt);
        simulate_short_wave!(canopy_rt, canOpt_rt, canRad_rt, sunRad_rt, soil);
        computeSIF_Fluxes!(arrayOfLeaves, canOpt_rt, canRad_rt, canopy_rt, soil, wl_set);
        push!(reflVIS, canRad_rt.alb_obs[28])
        push!(reflNIR, canRad_rt.alb_obs[52])
        push!(SIF_R , canRad_rt.SIF_obs[8])
        push!(SIF_FR, canRad_rt.SIF_obs[20])
    end
end
#----------------------------------------------------------------------------

A = reshape(reflNIR, ( 86,361));
B = reshape(reflVIS, ( 86,361));
SIFFER = reshape(SIF_R, ( 86,361));
SIFFER_FR = reshape(SIF_FR, ( 86,361));
#----------------------------------------------------------------------------

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()

# ## Setup
using Dojo

# ## Visualizer
vis = Visualizer()
open(vis)

################################################################################
## Astronaut
################################################################################
## humanoid
## multiple bodies
## initial linear velocities
## no gravity
## no spring and damper
## no control
################################################################################

function astronaut_simulation(vis::Visualizer;
        timestep=1.0e-2,
        tsim=1.0,
        gravity=0.0,
        spring=0.0,
        damper=0.0,
        seed=100,
        ϵ=1.0e-14,
        display=true)

    Random.seed!(seed)
    mech = get_mechanism(:humanoid,
        timestep=timestep,
        gravity=gravity,
        spring=spring,
        damper=damper,
        contact=false)
    initialize!(mech, :humanoid)

    ## Initialize bodies with random velocities
    for body in mech.bodies
        v = rand(3) .- 0.5
        v ./= norm(v)
        set_maximal_velocities!(body,
            v=v)
    end

    tcompute = @elapsed storage = simulate!(mech, tsim,
        record=true,
        opts=SolverOptions(rtol=ϵ, btol=ϵ))

    return mech, storage, tcompute
end

# ## Simulation
solver_tolerance = 1.0e-14
timestep = [0.05, 0.01, 0.005]
tsim = 10.0
storage = []
mech, _ = astronaut_simulation(vis,
    tsim=tsim,
    timestep=timestep[1],
    ϵ=solver_tolerance)

for t in timestep
    mech, data, tcompute = astronaut_simulation(vis,
        tsim=tsim,
        timestep=t,
        ϵ=solver_tolerance)
    push!(storage, data)
end

# ## Energy
color = [:magenta, :orange, :cyan];
me = [Dojo.mechanical_energy(mech, s)[2:end] for s in storage]
plt = plot(xlabel="time (s)", ylabel="energy drift");
for (i, e) in enumerate(me)
    tt = [j * timestep[i] for j = 1:length(e)]
    plt = plot!(tt, e .- e[1], color=color[i], width=2.0, label="h = $(timestep[i])")
end
display(plt);

# ## Momentum
m = [Dojo.momentum(mech, s)[1:end] for s in storage]

## linear
mlin = [[Vector(mi)[1:3] for mi in mt] for mt in m]
plt = plot(ylims=(-1.0e-10, 1.0e-10), xlabel="time (s)", ylabel="linear momentum drift");
for (i, mt) in enumerate(mlin)
    tt = [j * timestep[i] for j in 1:length(mt)]
    plt = plot!(tt, hcat([mi - mt[1] for mi in mt]...)', color=color[i], label=["h = $(timestep[i])" "" ""])
end
display(plt)

## angular
mang = [[Vector(mi)[4:6] for mi in mt] for mt in m]
plt = plot(ylims=(-1.0e-10, 1.0e-10), xlabel="time (s)", ylabel="angular momentum drift");
for (i, mt) in enumerate(mang)
    tt = [j * timestep[i] for j in 1:length(mt)]
    plt = plot!(tt, hcat([mi - mt[1] for mi in mt]...)', color=color[i], label=["h = $(timestep[i])" "" ""])
end
display(plt)


text:
using Test, AlphaVantage, Dates, DataFrames, CryptoDashApp, HTTP, JSON, Query

AlphaVantage.global_key!("AKTJ25ALEBBLH1QJ")

# Check if "data" folder exists, if not, create a new one
if isdir("data")
    @info "data folder exists, cleanup action will be performed!"
else
    mkdir("data")
    @info "New data folder has been created"
end

# Perform cleanup
CryptoDashApp.remove_old_files()

################# Test cases for accessing market data #################

@testset "Check if AV market data are accessible" begin

    for currency in ["ETH", "LTC", "LINK"]

        df_out_price, df_out_candle, df_out_vol = CryptoDashApp.get_price_data_single(currency)

        @test ~isempty(df_out_price)
        @test ~isempty(df_out_candle)
        @test ~isempty(df_out_vol)

    end

end

@testset "Check for exception handling while accessing AV market data" begin

    currency = "dummy"

    @test_logs (:info, "Could not fetch data, try again later!") match_mode=:any CryptoDashApp.get_price_data_single(currency)

end

################# Test cases for moving averages #################

@testset "Check if MA, MACD and signal are calculated" begin

    for currency in ["BTC", "DOT"]

        df_out_price, _ , _ = CryptoDashApp.get_price_data_single(currency)

        duration = 90
        window = 30

        Price_SMA, Price_WMA, Price_EMA = CryptoDashApp.moving_averages(df_out_price, duration, window)

        @test ~isempty(Price_SMA)
        @test ~isempty(Price_WMA)
        @test ~isempty(Price_WMA)

        df_out_price = df_out_price[end-duration+1-26-9+1:end, :]
        df_ema_all = CryptoDashApp.calculate_macd(df_out_price)

        @test ~isempty(df_ema_all)

    end

end

################# Test cases for CoinGecko API  #################

@testset "Check if CG developer and community data are accessible" begin

    for currency in ["btc", "eth", "ltc"]

        df_dev, df_comm = CryptoDashApp.get_dev_comm_data(currency)

        @test ~isempty(df_dev)
        @test ~isempty(df_comm)

    end

end

@testset "Check for exception handling while determining coin id" begin

    currency = "dummy"

    @test_logs (:info, "Could not find an id for the given currency") match_mode=:any CryptoDashApp.get_coin_id(currency)

end

@testset "Check if CG exchange volume data per currency are accessible" begin

    for currency in ["btc", "eth", "ltc"]

        num_exchanges = 10

        df_ex_vol = CryptoDashApp.get_exchange_vol_data(currency, num_exchanges)

        @test size(df_ex_vol)[1] == num_exchanges

    end

end

@testset "Check if CG overall exchange volume data are accessible" begin

    num_exchanges = 10

    for duration in [5, 10, 50, 75]
        df_ex_vol = CryptoDashApp.get_overall_vol_data(duration, num_exchanges)

        # Check for rows
        @test size(df_ex_vol)[1] == duration

        # Check for columns, total is num_exchanges + 1 due to "Time" column
        @test size(df_ex_vol)[2] == num_exchanges + 1
    end

end

#=@testset "Check if ratings data is accessible" begin

    utility_score, fcas_score, dev_score, mark_score, fcas_rating = CryptoDashApp.get_ratings_data("ETH")

    @test ~isempty(utility_score)
    @test ~isempty(fcas_score)
    @test ~isempty(dev_score)
    @test ~isempty(mark_score)
    @test ~isempty(fcas_rating)
end=#


text:
module CuArrays

using CUDAapi, CUDAdrv, CUDAnative

using GPUArrays

export CuArray, CuVector, CuMatrix, CuVecOrMat, cu
export CUBLAS, CUSPARSE, CUSOLVER, CUFFT, CURAND, CUDNN, CUTENSOR

import LinearAlgebra

using Adapt

using Libdl

using Requires

## source code includes

include("memory.jl")
include("array.jl")
include("subarray.jl")
include("permuteddimsarray.jl")
include("utils.jl")
include("indexing.jl")
include("broadcast.jl")
include("matmul.jl")
include("mapreduce.jl")
include("accumulate.jl")
include("linalg.jl")

include("gpuarray_interface.jl")

# many libraries need to be initialized per-device (per-context, really, but we assume users
# of CuArrays and/or CUDAnative only use a single context), so keep track of the active one.
const active_context = Ref{CuContext}()

include("blas/CUBLAS.jl")
include("sparse/CUSPARSE.jl")
include("solver/CUSOLVER.jl")
include("fft/CUFFT.jl")
include("rand/CURAND.jl")
include("dnn/CUDNN.jl")
include("tensor/CUTENSOR.jl")

include("nnlib.jl")

include("deprecated.jl")

## initialization

const __initialized__ = Ref(false)
functional() = __initialized__[]

export has_cudnn, has_cutensor
const libraries = Dict{String,Union{String,Nothing}}()
has_cudnn() = libraries["cudnn"] !== nothing && CUDNN.libcudnn !== nothing
has_cutensor() = libraries["cutensor"] !== nothing && CUTENSOR.libcutensor !== nothing

function __init__()
    silent = parse(Bool, get(ENV, "JULIA_CUDA_SILENT", "false"))
    verbose = parse(Bool, get(ENV, "JULIA_CUDA_VERBOSE", "false"))
    precompiling = ccall(:jl_generating_output, Cint, ()) != 0

    # if any dependent GPU package failed, expect it to have logged an error and bail out
    if !CUDAdrv.functional() || !CUDAnative.functional()
        verbose && @warn "CuArrays.jl did not initialize because CUDAdrv.jl or CUDAnative.jl failed to"
        return
    end

    # if any dependent GPU package failed, expect it to have logged an error and bail out
    CUDAdrv.functional() || return
    CUDAnative.functional() || return

    precompiling = ccall(:jl_generating_output, Cint, ()) != 0
    try
        # discover libraries
        toolkit = find_toolkit()
        for name in ("cublas", "cusparse", "cusolver", "cufft", "curand", "cudnn", "cutensor")
            mod = getfield(CuArrays, Symbol(uppercase(name)))
            lib = Symbol("lib$name")
            path = find_cuda_library(name, toolkit)
            libraries[name] = path
            if path !== nothing
                dir = dirname(path)
                if !(dir in Libdl.DL_LOAD_PATH)
                    push!(Libdl.DL_LOAD_PATH, dir)
                end
            end
        end

        # library dependencies
        CUBLAS.version()
        CUSPARSE.version()
        CUSOLVER.version()
        CUFFT.version()
        CURAND.version()
        # CUDNN and CUTENSOR are optional

        # library compatibility
        if has_cutensor()
            ver = CUTENSOR.version()
            if ver.major != 1 || ver.minor != 0
                error("CuArrays.jl only supports CUTENSOR 1.0")
            end
        end

        # package integrations
        @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" include("forwarddiff.jl")

        # update the active context when we switch devices
        callback = (::CuDevice, ctx::CuContext) -> begin
            active_context[] = ctx

            # wipe the active handles
            CUBLAS._handle[] = C_NULL
            CUBLAS._xt_handle[] = C_NULL
            CUSOLVER._dense_handle[] = C_NULL
            CUSOLVER._sparse_handle[] = C_NULL
            CUSPARSE._handle[] = C_NULL
            CURAND._generator[] = nothing
            CUDNN._handle[] = C_NULL
            CUTENSOR._handle[] = C_NULL
        end
        push!(CUDAnative.device!_listeners, callback)

        # a device might be active already
        existing_ctx = CUDAdrv.CuCurrentContext()
        if existing_ctx !== nothing
            active_context[] = existing_ctx
        end

        __init_memory__()
    catch ex
        # don't actually fail to keep the package loadable
        silent = parse(Bool, get(ENV, "JULIA_CUDA_SILENT", "false"))
        if !silent && !precompiling
            verbose = parse(Bool, get(ENV, "JULIA_CUDA_VERBOSE", "false"))
            if verbose
                @error "CuArrays.jl failed to initialize" exception=(ex, catch_backtrace())
            else
                @info "CuArrays.jl failed to initialize and will be unavailable (set JULIA_CUDA_SILENT or JULIA_CUDA_VERBOSE to silence or expand this message)"
            end
        end
    end
end

end # module


text:
using LinearAlgebra

mutable struct EarthSurrogate{X,Y,L,U,B,C,P,M,N,R,G,I,T} <: AbstractSurrogate
    x::X
    y::Y
    lb::L
    ub::U
    basis::B
    coeff::C
    penalty::P
    n_min_terms::M
    n_max_terms::N
    rel_res_error::R
    rel_GCV::G
    intercept::I
    maxiters::T
 end

 _hinge(x::Number,knot::Number) = max(0,x-knot)
 _hinge_mirror(x::Number,knot::Number) = max(0,knot-x)

 function _coeff_1d(x,y,basis)
     n = length(x)
     d = length(basis)
     X = zeros(eltype(x[1]),n,d)
     @inbounds for i = 1:n
         for j = 1:d
             X[i,j] = basis[j](x[i])
         end
     end
     return (X'*X)\(X'*y)
 end

 function _forward_pass_1d(x,y,n_max_terms,rel_res_error,maxiters)
     n = length(x)
     basis = Array{Function}(undef,0)
     current_sse = +Inf
     intercept = sum([y[i] for i =1:length(y)])/length(y)
     num_terms = 0
     pos_of_knot = 0
     iters = 0
     while num_terms < n_max_terms || iters > maxiters
         #Look for best addition:
         new_addition = false
         for i = 1:length(x)
            #Add or not add the knot var_i?
            var_i = x[i]
            new_basis = copy(basis)
            #select best new pair
            hinge1 = x-> _hinge(x,var_i)
            hinge2 = x-> _hinge_mirror(x,var_i)
            push!(new_basis,hinge1)
            push!(new_basis,hinge2)
            #find coefficients
            d = length(new_basis)
            X = zeros(eltype(x[1]),n,d)
            @inbounds for i = 1:n
                for j = 1:d
                    X[i,j] = new_basis[j](x[i])
                end
            end
            if (cond(X'*X) > 1e8)
                condition_number = false
                new_sse = +Inf
            else
                condition_number = true
                coeff = (X'*X)\(X'*y)
                new_sse = zero(y[1])
                d = length(new_basis)
                for i = 1:n
                    val_i = sum(coeff[j]*new_basis[j](x[i]) for j = 1:d) + intercept
                    new_sse = new_sse + (y[i]-val_i)^2
                end
            end
            #is the i-esim the best?
            if ( (new_sse < current_sse) && (abs(current_sse - new_sse) >= rel_res_error) && condition_number)
                #Add the hinge function to the basis
                pos_of_knot = i
                current_sse = new_sse
                new_addition = true
            end
        end
        iters = iters+1
        if new_addition
            best_hinge1 = z-> _hinge(z,x[pos_of_knot])
            best_hinge2 = z-> _hinge_mirror(z,x[pos_of_knot])
            push!(basis,best_hinge1)
            push!(basis,best_hinge2)
            num_terms = num_terms + 1
        else
            break
        end
     end
    if length(basis) == 0
        throw("Earth surrogate did not add any term, just the intercept. It is advised to double check the parameters.")
    end
    return basis
 end

 function _backward_pass_1d(x,y,n_min_terms,basis,penalty,rel_GCV)
     n = length(x)
     d = length(basis)
     intercept = sum([y[i] for i =1:length(y)])/length(y)
     coeff = _coeff_1d(x,y,basis)
     sse = zero(y[1])
     for i = 1:n
         val_i = sum(coeff[j]*basis[j](x[i]) for j = 1:d) + intercept
         sse = sse + (y[i]-val_i)^2
     end
     effect_num_params = d + penalty*(d-1)/2
     current_gcv = sse/(n*(1-effect_num_params/n)^2)
     num_terms = d
     while (num_terms > n_min_terms)
         #Basis-> select worst performing element-> eliminate it
         if num_terms <= 1
             break
         end
         found_new_to_eliminate = false
         for i = 1:n
             current_basis = copy(basis)
             #remove i-esim element from current basis
             deleteat!(current_basis,i)
             coef = _coeff_1d(x,y,current_basis)
             new_sse = zero(y[i])
             for a = 1:n
                 val_a = sum(coeff[j]*basis[j](x[a]) for j = 1:num_terms) + intercept
                 new_sse = new_sse + (y[a]-val_a)^2
             end
             effect_num_params = num_terms + penalty*(num_terms-1)/2
             i_gcv = new_sse/(n*(1-effect_num_params/n)^2)
             if i_gcv < current_gcv
                 basis_to_remove = i
                 new_gcv = i_gcv
                 found_new_to_eliminate = true
             end
         end
         if !found_new_to_eliminate
             break
         end
         if abs(current_gcv-new_gcv) < rel_GCV
             break
         else
             num_terms = num_terms-1
             deleteat!(basis,basis_to_remove)
         end
     end
     return basis
 end

function EarthSurrogate(x,y,lb::Number,ub::Number; penalty::Number = 2.0, n_min_terms::Int = 2, n_max_terms::Int = 10, rel_res_error::Number = 1e-2, rel_GCV::Number = 1e-2, maxiters = 100)
    intercept = sum([y[i] for i =1:length(y)])/length(y)
    basis_after_forward = _forward_pass_1d(x,y,n_max_terms,rel_res_error,maxiters)
    basis = _backward_pass_1d(x,y,n_min_terms,basis_after_forward,penalty,rel_GCV)
    coeff = _coeff_1d(x,y,basis)
    return EarthSurrogate(x,y,lb,ub,basis,coeff,penalty,n_min_terms,n_max_terms,rel_res_error,rel_GCV,intercept,maxiters)
end

function (earth::EarthSurrogate)(val::Number)
    return sum([earth.coeff[i]*earth.basis[i](val) for i = 1:length(earth.coeff)])+earth.intercept
end

 #ND
 #inside arr_hing I have functions like g(x) = x -> _hinge(x,5.0) or g(x) = one(x)
#_product_hinge(val,arr_hing) = prod([arr_hing[i](val[i]) for i = 1:length(val)])

function _coeff_nd(x,y,basis)
    n = length(x)
    base_len = length(basis)
    d = length(x[1])
    X = zeros(eltype(x[1]),n,base_len)
    @inbounds for a = 1:n
        for b = 1:base_len
            X[a,b] = prod([basis[b][c](x[a][c]) for c = 1:d])
        end
    end
    return (X'*X)\(X'*y)
end

function _forward_pass_nd(x,y,n_max_terms,rel_res_error,maxiters)
    n = length(x)
    basis = Array{Array{Function,1}}(undef,0) #ex of basis: push!(x,[x->1.0,x->1.0,x->_hinge(x,knot)]) so basis is of the form arr_hing and then I can pass my val
    current_sse = +Inf
    const_1 = x->1.0
    intercept = sum([y[i] for i =1:length(y)])/length(y)
    num_terms = 0
    d = length(x[1])
    best_hinge1 = best_hinge2 = [x->one(eltype(x[1])) for j = 1:d]
    new_addition = false
    iters = 0
    while num_terms <= n_max_terms || iters > maxiters
        current_basis = copy(basis)
        new_addition = false
        for i = 1:n
            for j = 1:d
                for k = 1:n
                    for l = 1:d
                        new_basis = copy(basis)
                        #model with interaction between x[i][j] and x[k][l]
                        new_hinge1 = [w != j ? x->one(eltype(x[1])) : z->_hinge(z,x[i][j]) for w = 1:d]
                        new_hinge2 = [w != l ? x->one(eltype(x[1])) : z->_hinge_mirror(z,x[k][l]) for w = 1:d]
                        push!(new_basis,new_hinge1)
                        push!(new_basis,new_hinge2)

                        #build the model, find sse and check if it's better than current best, if so save the params
                        bas_len = length(new_basis)
                        X = zeros(eltype(x[1]),n,bas_len)
                        @inbounds for a = 1:n
                            for b = 1:bas_len
                                X[a,b] = prod([new_basis[b][c](x[a][c]) for c = 1:d])
                            end
                        end
                        if (cond(X'*X) > 1e8)
                            condition_number = false
                            new_sse = +Inf
                        else
                            condition_number = true
                            coeff = (X'*X)\(X'*y)
                            new_sse = zero(y[1])
                            for a = 1:n
                                val_a = sum(coeff[b]*prod([new_basis[b][c](x[a][c]) for c = 1:d]) for b = 1:bas_len) + intercept
                                new_sse = new_sse + (y[a]-val_a)^2
                            end
                        end
                        #is the i-esim the best?
                        if ( (new_sse < current_sse) && (abs(current_sse - new_sse) >= rel_res_error) && condition_number)
                            #Add the hinge function to the basis
                            best_hinge1 = new_hinge1
                            best_hinge2 = new_hinge2
                            current_sse = new_sse
                            new_addition = true
                        end
                    end

                end
            end
        end
        iters = iters+1
        if new_addition
            push!(basis,best_hinge1)
            push!(basis,best_hinge2)
            num_terms = num_terms + 1
        else
            break
        end
    end
    if (length(basis) ==0)
        throw("Earth surrogate did not add any term, just the intercept. It is advised to double check the parameters.")
    end
    return basis
end

function _backward_pass_nd(x,y,n_min_terms,basis,penalty,rel_GCV)
    n = length(x)
    d = length(x[1])
    base_len = length(basis)
    intercept = sum([y[i] for i =1:length(y)])/length(y)
    coeff = _coeff_nd(x,y,basis)
    sse = zero(y[1])
    for a = 1:n
        val_a = sum(coeff[b]*prod([basis[b][c](x[a][c]) for c = 1:d]) for b = 1:base_len) + intercept
        sse = sse + (y[a]-val_a)^2
    end
    effect_num_params = base_len + penalty*(base_len-1)/2
    current_gcv = sse/(n*(1-effect_num_params/n)^2)
    num_terms = base_len
    new_gcv = +Inf
    basis_to_remove = 0
    while (num_terms > n_min_terms)
        #Basis-> select worst performing element-> eliminate it
        if num_terms <= 1
            break
        end
        found_new_to_eliminate = false
        for i = 1:num_terms
            current_basis = copy(basis)
            #remove i-esim element from current basis
            deleteat!(current_basis,i)
            coef = _coeff_nd(x,y,current_basis)
            new_sse = zero(y[i])
            current_base_len = num_terms-1
            for a = 1:n
                val_a = sum(coeff[b]*prod([current_basis[b][c](x[a][c]) for c = 1:d]) for b = 1:current_base_len) + intercept
                new_sse = new_sse + (y[a]-val_a)^2
            end
            curr_effect_num_params = current_base_len + penalty*(current_base_len - 1)/2
            i_gcv = new_sse/(n*(1-curr_effect_num_params/n)^2)
            if i_gcv < current_gcv
                basis_to_remove = i
                new_gcv = i_gcv
                found_new_to_eliminate = true
            end
        end
        if !found_new_to_eliminate
            break
        elseif abs(current_gcv-new_gcv) < rel_GCV
            break
        else
            num_terms = num_terms-1
            deleteat!(basis,basis_to_remove)
        end
    end
    return basis
end

function EarthSurrogate(x,y,lb,ub; penalty::Number = 2.0, n_min_terms::Int = 2, n_max_terms::Int = 10, rel_res_error::Number = 1e-2, rel_GCV::Number = 1e-2, maxiters = 100)
    intercept = sum([y[i] for i =1:length(y)])/length(y)
    basis_after_forward = _forward_pass_nd(x,y,n_max_terms,rel_res_error,maxiters)
    basis = _backward_pass_nd(x,y,n_min_terms,basis_after_forward,penalty,rel_GCV)
    coeff = _coeff_nd(x,y,basis)
    return EarthSurrogate(x,y,lb,ub,basis,coeff,penalty,n_min_terms,n_max_terms,rel_res_error,rel_GCV,intercept,maxiters)
end

function (earth::EarthSurrogate)(val)
    return sum([earth.coeff[i]*prod([earth.basis[i][j](val[j]) for j = 1:length(val)]) for i = 1:length(earth.coeff)]) + earth.intercept
end

function add_point!(earth::EarthSurrogate,x_new,y_new)
      if length(earth.x[1]) == 1
          #1D
          earth.x = vcat(earth.x,x_new)
          earth.y = vcat(earth.y,y_new)
          earth.intercept = sum([earth.y[i] for i =1:length(earth.y)])/length(earth.y)
          basis_after_forward = _forward_pass_1d(earth.x,earth.y,earth.n_max_terms,earth.rel_res_error,earth.maxiters)
          earth.basis = _backward_pass_1d(earth.x,earth.y,earth.n_min_terms,basis_after_forward,earth.penalty,earth.rel_GCV)
          earth.coeff = _coeff_1d(earth.x,earth.y,earth.basis)
          nothing
      else
          #ND
          earth.x = vcat(earth.x,x_new)
          earth.y = vcat(earth.y,y_new)
          earth.intercept = sum([earth.y[i] for i =1:length(earth.y)])/length(earth.y)
          basis_after_forward = _forward_pass_nd(earth.x,earth.y,earth.n_max_terms,earth.rel_res_error,earth.maxiters)
          earth.basis = _backward_pass_nd(earth.x,earth.y,earth.n_min_terms,basis_after_forward,earth.penalty,earth.rel_GCV)
          earth.coeff = _coeff_nd(earth.x,earth.y,earth.basis)
          nothing
      end

  end


text:
using Flux: onehotbatch, OneHotArray

function get_category_dict(X_cat::Array)
    category_dict = Dict{Float64, Vector{Float32}}()
    for variable in 1:size(X_cat, 1)
        category_dict[variable |> Float64] = unique(X_cat[variable, :])
    end
    return category_dict
end

"""
Build categorical feature embedding using
one-hot encoding
"""
function ohe_cat_features(X_cat, category_dict)

    # loop over columns and build tuple of ohe matricies
    X_cat_ohe =  Any[]
    row_ix = 1
    @inbounds for row in eachrow(X_cat)
        push!(X_cat_ohe, onehotbatch(row, category_dict[row_ix]))
        row_ix += 1
    end

    return X_cat_ohe
end


text:
# MLJ interface for Decision Trees/Random Forests models

import MLJModelInterface       # It seems that having done this in the top module is not enought
const MMI = MLJModelInterface  # We need to repoeat it here

export PerceptronClassifier, KernelPerceptronClassifier, PegasosClassifier

# ------------------------------------------------------------------------------
# Model Structure declarations..

mutable struct PerceptronClassifier <: MMI.Probabilistic
   initialθ::Vector{Float64}
   initialθ₀::Float64
   maxEpochs::Int64
   shuffle::Bool
   forceOrigin::Bool
   returnMeanHyperplane::Bool
   rng::AbstractRNG
end
PerceptronClassifier(;
  initialθ=Float64[],
  initialθ₀=0.0,
  maxEpochs=1000,
  shuffle=false,
  forceOrigin=false,
  returnMeanHyperplane=false,
  rng = Random.GLOBAL_RNG,
  ) = PerceptronClassifier(initialθ,initialθ₀,maxEpochs,shuffle,forceOrigin,returnMeanHyperplane,rng)

#kernelPerceptron(x, y; K=radialKernel, T=100, α=zeros(Int64,length(y)), nMsgs=0, shuffle=false)
 mutable struct KernelPerceptronClassifier <: MMI.Probabilistic
     K::Function
     maxEpochs::Int64
     initialα::Vector{Int64}
     shuffle::Bool
     rng::AbstractRNG
 end
 KernelPerceptronClassifier(;
    K=radialKernel,
    maxEpochs=100,
    initialα = Int64[],
    shuffle=false,
    rng = Random.GLOBAL_RNG,
    ) = KernelPerceptronClassifier(K,maxEpochs,initialα,shuffle,rng)

# pegasos(x, y; θ=zeros(size(x,2)),θ₀=0.0, λ=0.5,η= (t -> 1/sqrt(t)), T=1000, nMsgs=10, shuffle=false,forceOrigin=false,returnMeanHyperplane=false
mutable struct PegasosClassifier <: MMI.Probabilistic
   initialθ::Vector{Float64}
   initialθ₀::Float64
   λ::Float64
   η::Function
   maxEpochs::Int64
   shuffle::Bool
   forceOrigin::Bool
   returnMeanHyperplane::Bool
   rng::AbstractRNG
end
PegasosClassifier(;
  initialθ=Float64[],
  initialθ₀=0.0,
  λ = 0.5,
  η = (t -> 1/sqrt(t)),
  maxEpochs=1000,
  shuffle=false,
  forceOrigin=false,
  returnMeanHyperplane=false,
  rng = Random.GLOBAL_RNG,
  ) = PegasosClassifier(initialθ,initialθ₀,λ,η,maxEpochs,shuffle,forceOrigin,returnMeanHyperplane,rng)

# ------------------------------------------------------------------------------
# Fit functions...

function MMI.fit(model::PerceptronClassifier, verbosity, X, y)
 x = MMI.matrix(X)                     # convert table to matrix
 allClasses = levels(y)
 initialθ  = length(model.initialθ) == 0 ? zeros(size(x,2)) : model.initialθ
 fitresult = perceptron(x, y; θ=initialθ, θ₀=model.initialθ₀, T=model.maxEpochs, nMsgs=0, shuffle=model.shuffle, forceOrigin=model.forceOrigin, returnMeanHyperplane=model.returnMeanHyperplane,rng=model.rng)
 cache=nothing
 report=nothing
 return (fitresult,allClasses), cache, report
end

function MMI.fit(model::KernelPerceptronClassifier, verbosity, X, y)
 x          = MMI.matrix(X)                     # convert table to matrix
 allClasses = levels(y)
 initialα   = length(model.initialα) == 0 ? zeros(Int64,length(y)) : model.initialα
 fitresult  = kernelPerceptron(x, y; K=model.K, T=model.maxEpochs, α=initialα, nMsgs=0, shuffle=model.shuffle,rng=model.rng)
 cache      = nothing
 report     = nothing
 return (fitresult,allClasses), cache, report
end

function MMI.fit(model::PegasosClassifier, verbosity, X, y)
 x = MMI.matrix(X)                     # convert table to matrix
 allClasses = levels(y)
 initialθ  = length(model.initialθ) == 0 ? zeros(size(x,2)) : model.initialθ
 fitresult = pegasos(x, y; θ=initialθ,θ₀=model.initialθ₀, λ=model.λ,η=model.η, T=model.maxEpochs, nMsgs=0, shuffle=model.shuffle, forceOrigin=model.forceOrigin, returnMeanHyperplane=model.returnMeanHyperplane,rng=model.rng)
 cache=nothing
 report=nothing
 return (fitresult,allClasses), cache, report
end

# ------------------------------------------------------------------------------
# Predict functions....

function MMI.predict(model::Union{PerceptronClassifier,PegasosClassifier}, fitresult, Xnew)
    fittedModel      = fitresult[1]
    #classes          = CategoricalVector(fittedModel.classes)
    classes          = fittedModel.classes
    allClasses       = fitresult[2] # as classes do not includes classes unsees at training time
    nLevels          = length(allClasses)
    nRecords         = MMI.nrows(Xnew)
    modelPredictions = Perceptron.predict(MMI.matrix(Xnew), fittedModel.θ, fittedModel.θ₀, fittedModel.classes)
    predMatrix       = zeros(Float64,(nRecords,nLevels))
    # Transform the predictions from a vector of dictionaries to a matrix
    # where the rows are the PMF of each record
    for n in 1:nRecords
        for (c,cl) in enumerate(allClasses)
            predMatrix[n,c] = get(modelPredictions[n],cl,0.0)
        end
    end
    #predictions = [MMI.UnivariateFinite(classes, predMatrix[i,:])
    #               for i in 1:nRecords]
    predictions = MMI.UnivariateFinite(allClasses,predMatrix,pool=missing)
    return predictions
end

function MMI.predict(model::KernelPerceptronClassifier, fitresult, Xnew)
    fittedModel      = fitresult[1]
    #classes          = CategoricalVector(fittedModel.classes)
    classes          = fittedModel.classes
    allClasses       = fitresult[2] # as classes do not includes classes unsees at training time
    nLevels          = length(allClasses)
    nRecords         = MMI.nrows(Xnew)
    #ŷtrain = Perceptron.predict([10 10; 2.2 2.5],model.x,model.y,model.α, model.classes,K=model.K)
    modelPredictions = Perceptron.predict(MMI.matrix(Xnew), fittedModel.x, fittedModel.y, fittedModel.α, fittedModel.classes, K=fittedModel.K)
    predMatrix       = zeros(Float64,(nRecords,nLevels))
    # Transform the predictions from a vector of dictionaries to a matrix
    # where the rows are the PMF of each record
    for n in 1:nRecords
        for (c,cl) in enumerate(allClasses)
            predMatrix[n,c] = get(modelPredictions[n],cl,0.0)
        end
    end
    #predictions = [MMI.UnivariateFinite(classes, predMatrix[i,:])
    #              for i in 1:nRecords]
    #predictions = MMI.UnivariateFinite(classes, predMatrix)
    predictions = MMI.UnivariateFinite(allClasses,predMatrix,pool=missing)
    #predictions4 = MMI.UnivariateFinite(modelPredictions,pool=classes,ordered=false)
    #predictions = MMI.UnivariateFinite(modelPredictions,pool=fittedModel.classes)
    return predictions
end

# ------------------------------------------------------------------------------
# Model metadata for registration in MLJ...

MMI.metadata_model(PerceptronClassifier,
    input_scitype    = MMI.Table(MMI.Infinite),
    target_scitype   = AbstractVector{<: MMI.Finite},
    supports_weights = false,
    descr            = "The classical perceptron algorithm using one-vs-all for multiclass, from the Beta Machine Learning Toolkit (BetaML).",
	load_path        = "BetaML.Perceptron.PerceptronClassifier"
)

MMI.metadata_model(KernelPerceptronClassifier,
    input_scitype    = MMI.Table(MMI.Infinite),
    target_scitype   = AbstractVector{<: MMI.Finite},
    supports_weights = false,
    descr            = "The kernel perceptron algorithm using one-vs-one for multiclass, from the Beta Machine Learning Toolkit (BetaML).",
	load_path        = "BetaML.Perceptron.KernelPerceptronClassifier"
)

MMI.metadata_model(PegasosClassifier,
    input_scitype    = MMI.Table(MMI.Infinite),
    target_scitype   = AbstractVector{<: MMI.Finite},
    supports_weights = false,
    descr            = "The gradient-based linear \"pegasos\" classifier using one-vs-all for multiclass, from the Beta Machine Learning Toolkit (BetaML).",
	load_path        = "BetaML.Perceptron.PegasosClassifier"
)


text:
# script file to reproduce benchmarks in the analysis at
# TODO http ref
using SearchBenchmarks
using Test

path = "" # working directory including path delimiter at end
source = path*"somejuliasource.jl"

benchmark(200, 100000, source, [naivesearch,bloom_v0], path*"(1) TXT julia impl much better than naive search",false)
benchmark(200, 100000, source, [bloom_v0,bloom_v1,bloom_v2,bloom_v3], path*"(2) TXT julia bloom impl optimizations",false)
benchmark(200, 100000, source, [bloom_v0,bloom_v3,bloom_v4], path*"(3) TXT do bloom test first",false)
benchmark(2000, 100000, source, [naivesearch,bloom_v0,bloom_v3,bloom_v4], path*"(4) TXT degeneration on large patterns",false)
benchmark(2000, 100000, source, [naivesearch,bloom_v3,bloom_v4,bloom_best], path*"(5) TXT limited bits in bloom filter",false)
benchmark(200, 100000, source, [bloom_v0,bloom_v3,bloom_best,bmsearch], path*"(6) TXT boyer moore alternative",false)
benchmark(2000, 100000, source, [naivesearch,bloom_v0,bloom_v3,bloom_best,bmsearch], path*"(7) TXT boyer moore on large patterns",false)


text:
# Hypothesis Verification for 3D Object Recognition
# http://www.pointclouds.org/documentation/tutorials/global_hypothesis_verification.php#global-hypothesis-verification

using PCLCommon
using PCLIO
using PCLFilters
using PCLFeatures
using PCLKDTree
using PCLRecognition
using PCLRegistration
using PCLVisualization
using Cxx

import PCLCommon: CorrespondenceVal

model_pcd = Pkg.dir("PCLIO", "test", "data", "milk.pcd")
scene_pcd = Pkg.dir("PCLIO", "test", "data", "milk_cartoon_all_small_clorox.pcd")

T = PointXYZRGBA
NT = Normal
DT = SHOT352
RT = ReferenceFrame

model_ss = 0.02
scene_ss = 0.02
descr_rad = 0.02
rf_rad = 0.015
cg_size = 0.035
cg_thresh = 5.0
use_hough = true

# Hypothesis verification
hv_clutter_reg = 5.0
hv_inlier_th = 0.005
hv_occlusion_th = 0.01
hv_rad_clutter = 0.03
hv_regularizer = 3.0
hv_rad_normals = 0.05
hv_detect_clutter = false

# ICP
icp_max_iter = 5
icp_corr_distance = 0.005

model = PointCloud{T}(model_pcd)
scene = PointCloud{T}(scene_pcd)
model_keypoints = PointCloud{T}()
scene_keypoints = PointCloud{T}()
model_normals = PointCloud{NT}()
scene_normals = PointCloud{NT}()
model_descriptors = PointCloud{DT}()
scene_descriptors = PointCloud{DT}()

info("Normal estimation")
norm_est = NormalEstimationOMP{T,NT}()
setKSearch(norm_est, 10)
setInputCloud(norm_est, model)
compute(norm_est, model_normals)

setInputCloud(norm_est, scene)
compute(norm_est, scene_normals)

# Downsample
info("Uniform sampling")
uniform_sampling = UniformSampling{T}()
setInputCloud(uniform_sampling, model)
PCLFilters.setRadiusSearch(uniform_sampling, model_ss)
filter(uniform_sampling, model_keypoints)
println("Model total points: $(length(model)); Selectd Key points: $(length(model_keypoints))")

setInputCloud(uniform_sampling, scene)
PCLFilters.setRadiusSearch(uniform_sampling, scene_ss)
filter(uniform_sampling, scene_keypoints)
println("Scene total points: $(length(scene)); Selectd Key points: $(length(scene_keypoints))")

info("Associating a 3D descriptor to each model and scene keypoint")
descr_est = SHOTEstimationOMP{T,NT,DT}()
PCLFeatures.setRadiusSearch(descr_est, descr_rad)
setInputCloud(descr_est, model_keypoints)
setInputNormals(descr_est, model_normals)
setSearchSurface(descr_est, model)
compute(descr_est, model_descriptors)

setInputCloud(descr_est, scene_keypoints)
setInputNormals(descr_est, scene_normals)
setSearchSurface(descr_est, scene)
compute(descr_est, scene_descriptors)

info("Find Model-Scene Correspondences with KdTree")
match_search = KdTreeFLANN{DT}()
setInputCloud(match_search, model_descriptors)

model_scene_corrs = Correspondences()
for i in 0:length(scene_descriptors)-1
    neigh_indices = icxx"std::vector<int>(1);"
    neigh_sqr_dists = icxx"std::vector<float>(1);"
    p = icxx"$(scene_descriptors[i]).descriptor[0];"
    !isfinite(Cfloat(p)) && continue

    found_neighs = nearestKSearch(match_search, scene_descriptors[i],
        1, neigh_indices, neigh_sqr_dists)
    if found_neighs == 1 && (Float32(neigh_sqr_dists[0]) < Float32(0.25))
        corr = CorrespondenceVal(neigh_indices[0], i, neigh_sqr_dists[0])
        push!(model_scene_corrs, corr)
    end
end
@show length(model_scene_corrs)

info("Clustering")
# TODO: this is apperently cheating.. Re-allocations causes segfault, so reserve enough
# memory in advance
rototranslations = icxx"""
    auto v = std::vector<Eigen::Matrix4f, Eigen::aligned_allocator<Eigen::Matrix4f>>();
    v.reserve(10);
    return v;
"""
clustered_corrs = icxx"std::vector<pcl::Correspondences>();"

clusterer = use_hough ? Hough3DGrouping{T,T,RT,RT}() :
    GeometricConsistencyGrouping{T,T}()

if use_hough
    model_rf = PointCloud{RT}()
    scene_rf = PointCloud{RT}()
    rf_est = BOARDLocalReferenceFrameEstimation{T,NT,RT}()
    setFindHoles(rf_est, true)
    PCLFeatures.setRadiusSearch(rf_est, rf_rad)

    setInputCloud(rf_est, model_keypoints)
    setInputNormals(rf_est, model_normals)
    setSearchSurface(rf_est, model)
    compute(rf_est, model_rf)

    setInputCloud(rf_est, scene_keypoints)
    setInputNormals(rf_est, scene_normals)
    setSearchSurface(rf_est, scene)
    compute(rf_est, scene_rf)

    setHoughBinSize(clusterer, cg_size)
    setHoughThreshold(clusterer, cg_thresh)
    setUseInterpolation(clusterer, true)
    setUseDistanceWeight(clusterer, false)
else
    setGCSize(clusterer, cg_size)
    setGCThreshold(clusterer, cg_thresh)
end

setInputCloud(clusterer, model_keypoints)
use_hough && setInputRf(clusterer, model_rf)
setSceneCloud(clusterer, scene_keypoints)
use_hough && setSceneRf(clusterer, scene_rf)
setModelSceneCorrespondences(clusterer, model_scene_corrs)
recognize(clusterer, rototranslations, clustered_corrs)

@show length(rototranslations)

# error("")
# instances = CxxStd.StdVector{cxxt"pcl::PointCloud<$T>::ConstPtr"}()
instances = icxx"std::vector<pcl::PointCloud<$T>::ConstPtr>();"
for i in 0:length(rototranslations)-1
    rotated_model = PointCloud{T}()
    # transformPointCloud(model, rotated_model, rototranslations[i])
    @show i
    hoge = icxx"$(rototranslations)[$i];"
    transformPointCloud(model, rotated_model, hoge)
    push!(instances, rotated_model.handle)
end

registered_instances = icxx"std::vector<pcl::PointCloud<$T>::ConstPtr>();"

info("ICP")
for i in 0:length(rototranslations)-1
    icp = IterativeClosestPoint{T,T}()
    setMaximumIterations(icp, icp_max_iter)
    setMaxCorrespondenceDistance(icp, icp_corr_distance)
    setInputTarget(icp, scene)
    setInputSource(icp, instances[i])
    registered = PointCloud{T}()
    align(icp, registered)

    print("Instance $i ")
    if hasConverged(icp)
        println("Aligned!")
    else
        println("Not aligned!")
    end
    push!(registered_instances, registered.handle)
end

hypothesis_mask = icxx"std::vector<bool>();"

GoHv = GlobalHypothesesVerification{T,T}()
setSceneCloud(GoHv, scene)
addModels(GoHv, registered_instances, true)
setInlierThreshold(GoHv, hv_inlier_th)
setRegularizer(GoHv, hv_regularizer)
setRadiusClutter(GoHv, hv_rad_clutter)
setClutterRegularizer(GoHv, hv_clutter_reg)
setDetectClutter(GoHv, hv_detect_clutter)
setRadiusNormals(GoHv, hv_rad_normals)

verify(GoHv)
getMask(GoHv, hypothesis_mask)
for i in 0:length(hypothesis_mask)-1
    if icxx"$hypothesis_mask[$i] == true;"
        println("Instance $i is good!")
    else
        println("Instance $i is bad!")
    end
end

if isdefined(:vis) && vis
    info("Prepare PCL visualizer...")
    viewer = PCLVisualizer("pcl visualizer")
    color_handler = PointCloudColorHandlerRGBField(scene)
    addPointCloud(viewer, scene, color_handler, id="scene")

    for i in 0:length(rototranslations)-1
        rotated_model = PointCloud{T}()
        transformPointCloud(model, rotated_model, rototranslations[i])

        r,g,b = icxx"$hypothesis_mask[$i]==true;" ? (255,0,0) : (0,255,0)
        color_handler = PointCloudColorHandlerCustom(scene, r, g, b)
        addPointCloud(viewer, rotated_model, color_handler, id="instance $i")
    end
    spin(viewer)
    close(viewer)
end


text:
module Model

using Gen
using ..MangroveBase: ReviewInfo, Rating, RATINGS

@gen function subset(items::Vector{String})::Vector
  filter(item -> @trace(bernoulli(0.5), item), items)
end

@dist rating() = uniform_discrete(1, RATINGS)

@gen function mangrove_model(subjects::Set{String}, reviewers::Set{String})
  for sub in subjects
    @trace(rating(), :quality => sub)
    for kid in @trace(subset(collect(reviewers)), :selected => sub)
      @trace(rating(), ReviewInfo(sub, kid) => :rating)
    end
  end
end

end

max_stars_repo_name: Open-Reviews/Mangrove

max_issues_repo_name: Open-Reviews/Mangrove

max_forks_repo_name: Open-Reviews/Mangrove