text:
using DataFrames

"""
Write global parameter file.
"""
function write_global_param(path_sim, path_forcing, startyear, endyear, timestep, output_force, full_energy, output_binary)

    # Set output binary true if output force is true

    output_binary = output_force == "TRUE" ? "TRUE" : output_binary

    # Standard parameter part

    str_standard = """
    #######################################################################
    # VIC Model Parameters - 4.1.x
    #######################################################################
    # Id: global.param.sample,v 5.7.2.27 2012/02/07 19:03:02 vicadmin Exp
    #######################################################################
    # Simulation Parameters
    #######################################################################
    NLAYER		3	# number of soil layers
    NODES		3	# number of soil thermal nodes
    TIME_STEP 	$timestep	# model time step in hours (set to 24 if FULL_ENERGY = FALSE, set to < 24 if FULL_ENERGY = TRUE)
    SNOW_STEP	$timestep	# time step in hours for which to solve the snow model (should = TIME_STEP if TIME_STEP < 24)
    STARTYEAR	$startyear	# year model simulation starts
    STARTMONTH	01	# month model simulation starts
    STARTDAY	01 	# day model simulation starts
    STARTHOUR	00	# hour model simulation starts
    ENDYEAR 	$endyear	# year model simulation ends
    ENDMONTH	12	# month model simulation ends
    ENDDAY		31	# day model simulation ends
    FULL_ENERGY 	$full_energy	# TRUE = calculate full energy balance; FALSE = compute water balance only
    FROZEN_SOIL	FALSE	# TRUE = calculate frozen soils

    #######################################################################
    # Soil Temperature Parameters
    # VIC will choose appropriate values for QUICK_FLUX and IMPLICIT depending on values of FULL_ENERGY and FROZEN_SOIL; the user should only need to override VICs choices in special cases.
    # The other options in this section are only applicable when FROZEN_SOIL is TRUE and their values depend on the appliqqcation.
    #######################################################################
    #QUICK_FLUX	FALSE	# TRUE = use simplified ground heat flux method of Liang et al (1999); FALSE = use finite element method of Cherkauer et al (1999)
    #IMPLICIT	TRUE	# TRUE = use implicit solution for soil heat flux equation of Cherkauer et al (1999), otherwise uses original explicit solution.
    #QUICK_SOLVE	FALSE	# TRUE = Use Liang et al., 1999 formulation for iteration, but explicit finite difference method for final step.
    #NO_FLUX		FALSE	# TRUE = use no flux lower boundary for ground heat flux computation; FALSE = use constant flux lower boundary condition.  If NO_FLUX = TRUE, QUICK_FLUX MUST = FALSE
    #EXP_TRANS	FALSE	# TRUE = exponentially distributes the thermal nodes in the Cherkauer et al. (1999) finite difference algorithm, otherwise uses linear distribution
    #GRND_FLUX_TYPE	GF_410	# Options for ground flux:
    #			# GF_406 = use (flawed) formulas for ground flux, deltaH, and fusion from VIC 4.0.6 and earlier;
    #			# GF_410 = use formulas from VIC 4.1.0 (ground flux, deltaH, and fusion are correct; deltaH and fusion ignore surf_atten);
    #			# Default = GF_410
    #TFALLBACK	TRUE	# TRUE = when temperature iteration fails to converge, use previous time steps T value

    #######################################################################
    # Precip (Rain and Snow) Parameters
    # Generally these default values do not need to be overridden
    #######################################################################
    #SNOW_ALBEDO	USACE	# USACE = use traditional VIC algorithm based on US Army Corps of Engineers empirical snow albedo decay curves, using hard-coded dates for transitions from snow accumulation to melting; SUN1999 = use algorithm of Sun et al 1999, in which albedo decay depends on snow cold content (more appropriate for simulations outside the US).
    #SNOW_DENSITY	DENS_BRAS	# DENS_BRAS = use traditional VIC algorithm taken from Bras, 1990; DENS_SNTHRM = use algorithm taken from SNTHRM model.
    #BLOWING		FALSE	# TRUE = compute evaporative fluxes due to blowing snow
    #COMPUTE_TREELINE	FALSE	# Can be either FALSE or the id number of an understory veg class; FALSE = turn treeline computation off; VEG_CLASS_ID = replace any overstory veg types with the this understory veg type in all snow bands for which the average July Temperature <= 10 C (e.g. COMPUTE_TREELINE 10 replaces any overstory veg cover with class 10)
    #DIST_PRCP	FALSE	# TRUE = use distributed precipitation
    #PREC_EXPT	0.6	# exponent for use in distributed precipitation eqn (only used if DIST_PRCP is TRUE)
    #CORRPREC	FALSE	# TRUE = correct precipitation for gauge underqqcatch
    #MAX_SNOW_TEMP	0.5	# maximum temperature (C) at which snow can fall
    #MIN_RAIN_TEMP	-0.5	# minimum temperature (C) at which rain can fall

    #######################################################################
    # Turbulent Flux Parameters
    # Generally these default values do not need to be overridden
    #######################################################################
    #MIN_WIND_SPEED	0.1	# minimum allowable wind speed (m/s)
    #AERO_RESIST_CANSNOW	AR_406_FULL	# Options for aerodynamic resistance in snow-filled canopy:
    #			# AR_406 	= multiply by 10 for latent heat but do NOT multiply by 10 for sensible heat and do NOT apply stability correction (as in VIC 4.0.6); when no snow in canopy, use surface aero_resist for ET.
    #			# AR_406_LS 	= multiply by 10 for latent heat AND sensible heat and do NOT apply stability correction; when no snow in canopy, use surface aero_resist for ET.
    #			# AR_406_FULL 	= multiply by 10 for latent heat AND sensible heat and do NOT apply stability correction; additionally, always use overstory aero_resist for ET (as in 4.1.0).
    #			# AR_410 	= apply stability correction but do NOT multiply by 10 (as in VIC 4.1.0); additionally, always use overstory aero_resist for ET (as in 4.1.0).
    #			# Default 	= AR_406_FULL

    #######################################################################
    # Meteorological Forcing Disaggregation Parameters
    # Generally these default values do not need to be overridden
    #######################################################################
    #PLAPSE		TRUE	# This controls how VIC computes air pressure when air pressure is not supplied as an input forcing: TRUE = set air pressure to sea level pressure, lapsed to grid cell average elevation; FALSE = set air pressure to constant 95.5 kPa (as in all versions of VIC pre-4.1.1)
    #SW_PREC_THRESH		0	# Minimum daily precip [mm] that can cause dimming of incoming shortwave; default = 0.
    MTCLIM_SWE_CORR	FALSE    # This controls VICs estimates of incoming shortwave in the presence of snow; TRUE = adjust incoming shortwave for snow albedo effect; FALSE = do not adjust shortwave; default = TRUE
    #VP_ITER		VP_ITER_ALWAYS	# This controls VICs iteration between estimates of shortwave and vapor pressure:
    #			# VP_ITER_NEVER = never iterate; make estimates separately
    #			# VP_ITER_ALWAYS = always iterate once
    #			# VP_ITER_ANNUAL = iterate once for arid climates based on annual Precip/PET ratio
    #			# VP_ITER_CONVERGE = iterate until shortwave and vp stabilize
    #			# default = VP_ITER_ALWAYS
    #VP_INTERP	TRUE	# This controls sub-daily humidity estimates; TRUE = interpolate daily VP estimates linearly between sunrise of one day to the next; FALSE = hold VP constant for entire day
    #LW_TYPE		LW_TVA	# This controls the algorithm used to estimate clear-sky longwave radiation:
    #			# LW_TVA = Tennessee Valley Authority algorithm (1972)
    #			# other options listed in vicNl_def.h
    #			# default = LW_TVA (this was also the traditional VIC algorithm)
    #LW_CLOUD	LW_CLOUD_DEARDORFF	# This controls the algorithm used to estimate the influence of clouds on total longwave:
    #			# LW_CLOUD_BRAS = method from Bras textbook (this was the traditional VIC algorithm)
    #			# LW_CLOUD_DEARDORFF = method of Deardorff (1978)
    #			# default = LW_CLOUD_DEARDORFF
    OUTPUT_FORCE	$output_force

    #######################################################################
    # Miscellaneous Simulation Parameters
    # Generally these default values do not need to be overridden
    #######################################################################
    #CONTINUEONERROR	TRUE	# TRUE = if simulation aborts on one grid cell, continue to next grid cell

    #######################################################################
    # State Files and Parameters
    #######################################################################
    #INIT_STATE	/home/jmg/VIC-master/12.70/params/vic/state/STEHE.state_19481231	# Initial state path/file
    #STATENAME	/home/jmg/VIC-master/12.70/state/STEHE.state	# Output state file path/prefix.  The date (STATEYEAR,STATEMONTH,STATEDAY) will be appended to the prefix automatically in the format yyyymmdd.
    #STATEYEAR	2000	# year to save model state
    #STATEMONTH	12	# month to save model state
    #STATEDAY	31	# day to save model state
    #BINARY_STATE_FILE       FALSE	# TRUE if state file should be binary format; FALSE if ascii

    #######################################################################
    # Forcing Files and Parameters
    #
    #       All FORCING filenames are actually the pathname, and prefix
    #               for gridded data types: ex. DATA/forcing_
    #               Latitude and longitude index suffix is added by VIC
    #
    #	There must be 1 FORCE_TYPE entry for each variable (column) in the forcing file
    #
    #	If FORCE_TYPE is BINARY, each FORCE_TYPE must be followed by:
    #			SIGNED/UNSIGNED	SCALE_FACTOR
    #		For example (BINARY):
    #			FORCE_TYPE	PREC	UNSIGNED	40
    #		or (ASCII):
    #			FORCE_TYPE	PREC
    #######################################################################

    FORCING1	$(path_forcing)/forcing/data_
    GRID_DECIMAL	5
    N_TYPES		4
    FORCE_TYPE	PREC      SIGNED 100
    FORCE_TYPE	TMIN      SIGNED   100
    FORCE_TYPE	TMAX      SIGNED   100
    FORCE_TYPE	WIND      SIGNED  100
    FORCE_FORMAT	BINARY
    FORCE_ENDIAN    LITTLE
    FORCE_DT	24
    FORCEYEAR	$startyear
    FORCEMONTH	1
    FORCEDAY	1
    FORCEHOUR	01
    WIND_H          10.0
    MEASURE_H       2.0
    ALMA_INPUT	FALSE

    #######################################################################
    # Land Surface Files and Parameters
    #######################################################################
    SOIL            $(path_sim)/params/soil_param	# Soil parameter path/file
    ARC_SOIL        FALSE   # TRUE = read soil parameters from ARC/INFO ASCII grids
    #SOIL_DIR        (soil param directory)   # Directory containing ARC/INFO ASCII grids of soil parameters - only valid if ARC_SOIL is TRUE
    BASEFLOW	NIJSSEN2001	# ARNO = columns 5-8 are the standard VIC baseflow parameters; NIJSSEN2001 = columns 5-8 of soil file are baseflow parameters from Nijssen et al (2001)
    JULY_TAVG_SUPPLIED	FALSE	# TRUE = final column of the soil parameter file will contain average July air temperature, for computing treeline; this will be ignored if COMPUTE_TREELINE is FALSE; FALSE = compute the treeline based on the average July air temperature of the forcings over the simulation period
    ORGANIC_FRACT	FALSE	# TRUE = simulate organic soils; soil param file contains 3*Nlayer extra columns, listing for each layer the organic fraction, and the bulk density and soil particle density of the organic matter in the soil layer; FALSE = soil param file does not contain any information about organic soil, and organic fraction should be assumed to be 0
    VEGLIB	        $(path_sim)/params/veglib_param	# Veg library path/file
    VEGPARAM        $(path_sim)/params/veg_param	# Veg parameter path/file
    ROOT_ZONES      2	# Number of root zones (must match format of veg param file)
    VEGPARAM_LAI 	TRUE    # TRUE = veg param file contains LAI information; FALSE = veg param file does NOT contain LAI information
    LAI_SRC 	LAI_FROM_VEGLIB    # LAI_FROM_VEGPARAM = read LAI from veg param file; LAI_FROM_VEGLIB = read LAI from veg library file
    SNOW_BAND	1 $(path_sim)/params/elev_band  # Number of snow bands; if number of snow bands > 1, you must insert the snow band path/file after the number of bands (e.g. SNOW_BAND 5 my_path/my_snow_band_file)

    #######################################################################
    # Lake Simulation Parameters
    # These need to be un-commented and set to correct values only when running lake model (LAKES is not FALSE)
    #######################################################################
    #LAKES		(put lake parameter path/file here)	# Lake parameter path/file
    #LAKE_PROFILE	FALSE	# TRUE = User-specified depth-area parameters in lake parameter file; FALSE = VIC computes a parabolic depth-area profile
    #EQUAL_AREA	FALSE	# TRUE = grid cells are from an equal-area projection; FALSE = grid cells are on a regular lat-lon grid
    #RESOLUTION	0.125	# Grid cell resolution (degrees if EQUAL_AREA is FALSE, km^2 if EQUAL_AREA is TRUE); ignored if LAKES is FALSE

    #######################################################################
    # Output Files and Parameters
    #######################################################################
    RESULT_DIR      $(path_sim)/results/	# Results directory path
    OUT_STEP        24       # Output interval (hours); if 0, OUT_STEP = TIME_STEP
    SKIPYEAR 	0	# Number of years of output to omit from the output files
    COMPRESS	FALSE	# TRUE = compress input and output files when done
    BINARY_OUTPUT $output_binary	# TRUE = binary output files
    ALMA_OUTPUT	FALSE	# TRUE = ALMA-format output files; FALSE = standard VIC units
    MOISTFRACT 	FALSE	# TRUE = output soil moisture as volumetric fraction; FALSE = standard VIC units
    PRT_HEADER	FALSE   # TRUE = insert a header at the beginning of each output file; FALSE = no header
    PRT_SNOW_BAND   FALSE   # TRUE = write a snowband output file, containing band-specific values of snow variables; NOTE: this is ignored if N_OUTFILES is specified below.

    #######################################################################
    #
    # Output File Contents
    #
    # As of VIC 4.0.6 and 4.1.0, you can specify your output file names and
    # contents # in the global param file (see the README.txt file for more
    # information).
    #
    # If you do not specify file names and contents in the global param
    # file, VIC will produce the same set of output files that it has
    # produced in earlier versions, namely fluxes and snow files, plus
    # fdepth files if FROZEN_SOIL is TRUE and snowband files if
    # PRT_SNOW_BAND is TRUE.  These files will have the same contents and
    # format as in earlier versions.
    #
    # The OPTIMIZE and LDAS_OUTPUT options have been removed.  These
    # output configurations can be selected with the proper set of
    # instructions in the global param file.  (see the output.*.template
    # files included in this distribution for more information.)
    #
    # If you do specify the file names and contents in the global param file,
    # PRT_SNOW_BAND will have no effect.
    #
    # Format:
    #
    #   N_OUTFILES    <n_outfiles>
    #
    #   OUTFILE       <prefix>        <nvars>
    #   OUTVAR        <varname>       [<format>        <type>  <multiplier>]
    #   OUTVAR        <varname>       [<format>        <type>  <multiplier>]
    #   OUTVAR        <varname>       [<format>        <type>  <multiplier>]
    #
    #   OUTFILE       <prefix>        <nvars>
    #   OUTVAR        <varname>       [<format>        <type>  <multiplier>]
    #   OUTVAR        <varname>       [<format>        <type>  <multiplier>]
    #   OUTVAR        <varname>       [<format>        <type>  <multiplier>]
    #
    #
    # where
    #   <n_outfiles> = number of output files
    #   <prefix>     = name of the output file, NOT including latitude
    #                  and longitude
    #   <nvars>      = number of variables in the output file
    #   <varname>    = name of the variable (this must be one of the
    #                  output variable names listed in vicNl_def.h.)
    #   <format>     = (for ascii output files) fprintf format string,
    #                  e.g.
    #                    %.4f = floating point with 4 decimal places
    #                    %.7e = scientific notation w/ 7 decimal places
    #                    *    = use the default format for this variable
    #
    #   <format>, <type>, and <multiplier> are optional.  For a given
    #   variable, you can specify either NONE of these, or ALL of
    #   these.  If these are omitted, the default values will be used.
    #
    #   <type>       = (for binary output files) data type code.
    #                  Must be one of:
    #                    OUT_TYPE_DOUBLE = double-precision floating point
    #                    OUT_TYPE_FLOAT  = single-precision floating point
    #                    OUT_TYPE_INT    = integer
    #                    OUT_TYPE_USINT  = unsigned short integer
    #                    OUT_TYPE_SINT   = short integer
    #                    OUT_TYPE_CHAR   = char
    #                    *               = use the default type
    #   <multiplier> = (for binary output files) factor to multiply
    #                  the data by before writing, to increase precision.
    #                    *    = use the default multiplier for this variable
    #
    #######################################################################

    """

    # Write mtclim binary output

    str_mtclim = """

    N_OUTFILES    1

    OUTFILE       metdata        8

    OUTVAR		OUT_PREC	* OUT_TYPE_USINT	40
    OUTVAR		OUT_AIR_TEMP	* OUT_TYPE_SINT		100
    OUTVAR		OUT_SHORTWAVE	* OUT_TYPE_USINT	50
    OUTVAR		OUT_LONGWAVE	* OUT_TYPE_USINT	80
    OUTVAR		OUT_PRESSURE	* OUT_TYPE_USINT	100
    OUTVAR		OUT_QAIR	* OUT_TYPE_USINT	100000
    OUTVAR		OUT_VP		* OUT_TYPE_USINT	100
    OUTVAR		OUT_WIND	* OUT_TYPE_USINT	100

    """

    # Write results as binary output

    str_results = """

    N_OUTFILES    1

    OUTFILE       results        28

    OUTVAR      OUT_PREC            * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_AIR_TEMP        * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_RAINF           * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_SNOWF           * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_WIND            * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_RUNOFF          * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_BASEFLOW        * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_INFLOW          * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_EVAP            * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_SUB_SNOW        * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_TRANSP_VEG      * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_SUB_CANOP       * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_ZWT_LUMPED      * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_SOIL_WET        * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_SWE             * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_SNOW_COVER      * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_SNOW_CANOPY     * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_WDEW            * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_LATENT          * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_SENSIBLE        * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_NET_LONG        * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_NET_SHORT       * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_PET_H2OSURF     * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_PET_SHORT       * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_PET_NATVEG      * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_SOIL_LIQ        * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_EVAP_CANOP      * OUT_TYPE_DOUBLE		*
    OUTVAR      OUT_EVAP_BARE       * OUT_TYPE_DOUBLE		*

    """

    # Concatenate standard parameters with mtclim output

    if output_force == "FALSE" && output_binary == "FALSE"
        str = str_standard
    elseif output_force == "FALSE" && output_binary == "TRUE"
        str = str_standard * str_results
    elseif output_force == "TRUE"
        str = str_standard * str_mtclim
    end

    # Write the global parameters to file

    fid = open(joinpath(path_sim, "params/global_param"), "w")

    write(fid, str)

    close(fid)

    return nothing

end

"""
Read soil parameter file.
"""
function read_soil_params(path)

    df = CSV.File(joinpath(path, "params/soil_param"), delim = ' ', header = false) |> DataFrame

    colnames = [:run_cell, :gridcel, :lat, :lon, :infilt, :D1, :D2, :D3, :D4,
    :expt1, :expt2, :expt3, :Ksat1, :Ksat2, :Ksat3, :phi_s1, :phi_s2, :phi_s3,
    :init_moist1, :init_moist2, :init_moist3, :elev, :depth1, :depth2, :depth3,
    :avg_T, :dp, :bubble1, :bubble2, :bubble3, :quartz1, :quartz2, :quartz3,
    :bulk_density1, :bulk_density2, :bulk_density3, :soil_density1, :soil_density2,
    :soil_density3, :off_gmt, :Wcr_FRACT1, :Wcr_FRACT2, :Wcr_FRACT3, :Wpwp_FRACT1,
    :Wpwp_FRACT2, :Wpwp_FRACT3, :rough, :snow_rough, :annual_prec, :resid_moist1,
    :resid_moist2, :resid_moist3, :fs_active]

    df = df[1:length(colnames)]

    names!(df, colnames)

    return df

end

"""
Write soil parameter file.
"""
function write_soil_params(path_sim, df)

    # Write parameters to file

    fid = open(joinpath(path_sim, "params/soil_param"), "w")

    for irow in 1:size(df,1)

        if df[:Wpwp_FRACT1][irow] < df[:resid_moist1][irow] / (1.0 - df[:bulk_density1][irow]/df[:soil_density1][irow])
            df[:Wpwp_FRACT1][irow] = 1.01 * df[:resid_moist1][irow] / (1.0 - df[:bulk_density1][irow]/df[:soil_density1][irow])
        end

        if df[:Wpwp_FRACT2][irow] < df[:resid_moist2][irow] / (1.0 - df[:bulk_density2][irow]/df[:soil_density2][irow])
            df[:Wpwp_FRACT2][irow] = 1.01 * df[:resid_moist2][irow] / (1.0 - df[:bulk_density2][irow]/df[:soil_density2][irow])
        end

        if df[:Wpwp_FRACT3][irow] < df[:resid_moist3][irow] / (1.0 - df[:bulk_density3][irow]/df[:soil_density3][irow])
            df[:Wpwp_FRACT3][irow] = 1.01 * df[:resid_moist3][irow] / (1.0 - df[:bulk_density3][irow]/df[:soil_density3][irow])
        end

        write(fid, @sprintf("%.0f ", df[:run_cell][irow]))
        write(fid, @sprintf("%.0f ", df[:gridcel][irow]))
        write(fid, @sprintf("%.5f ", df[:lat][irow]))
        write(fid, @sprintf("%.5f ", df[:lon][irow]))
        write(fid, @sprintf("%.4f ", df[:infilt][irow]))
        write(fid, @sprintf("%.4f ", df[:D1][irow]))
        write(fid, @sprintf("%.4f ", df[:D2][irow]))
        write(fid, @sprintf("%.4f ", df[:D3][irow]))
        write(fid, @sprintf("%.4f ", df[:D4][irow]))
        write(fid, @sprintf("%.4f ", df[:expt1][irow]))
        write(fid, @sprintf("%.4f ", df[:expt2][irow]))
        write(fid, @sprintf("%.4f ", df[:expt3][irow]))
        write(fid, @sprintf("%.4f ", df[:Ksat1][irow]))
        write(fid, @sprintf("%.4f ", df[:Ksat2][irow]))
        write(fid, @sprintf("%.4f ", df[:Ksat3][irow]))
        write(fid, @sprintf("%.4f ", df[:phi_s1][irow]))
        write(fid, @sprintf("%.4f ", df[:phi_s2][irow]))
        write(fid, @sprintf("%.4f ", df[:phi_s3][irow]))
        write(fid, @sprintf("%.4f ", df[:init_moist1][irow]))
        write(fid, @sprintf("%.4f ", df[:init_moist2][irow]))
        write(fid, @sprintf("%.4f ", df[:init_moist3][irow]))
        write(fid, @sprintf("%.4f ", df[:elev][irow]))
        write(fid, @sprintf("%.4f ", df[:depth1][irow]))
        write(fid, @sprintf("%.4f ", df[:depth2][irow]))
        write(fid, @sprintf("%.4f ", df[:depth3][irow]))
        write(fid, @sprintf("%.4f ", df[:avg_T][irow]))
        write(fid, @sprintf("%.4f ", df[:dp][irow]))
        write(fid, @sprintf("%.4f ", df[:bubble1][irow]))
        write(fid, @sprintf("%.4f ", df[:bubble2][irow]))
        write(fid, @sprintf("%.4f ", df[:bubble3][irow]))
        write(fid, @sprintf("%.4f ", df[:quartz1][irow]))
        write(fid, @sprintf("%.4f ", df[:quartz2][irow]))
        write(fid, @sprintf("%.4f ", df[:quartz3][irow]))
        write(fid, @sprintf("%.4f ", df[:bulk_density1][irow]))
        write(fid, @sprintf("%.4f ", df[:bulk_density2][irow]))
        write(fid, @sprintf("%.4f ", df[:bulk_density3][irow]))
        write(fid, @sprintf("%.4f ", df[:soil_density1][irow]))
        write(fid, @sprintf("%.4f ", df[:soil_density2][irow]))
        write(fid, @sprintf("%.4f ", df[:soil_density3][irow]))
        write(fid, @sprintf("%.4f ", df[:off_gmt][irow]))
        write(fid, @sprintf("%.4f ", df[:Wcr_FRACT1][irow]))
        write(fid, @sprintf("%.4f ", df[:Wcr_FRACT2][irow]))
        write(fid, @sprintf("%.4f ", df[:Wcr_FRACT3][irow]))
        write(fid, @sprintf("%.4f ", df[:Wpwp_FRACT1][irow]))
        write(fid, @sprintf("%.4f ", df[:Wpwp_FRACT2][irow]))
        write(fid, @sprintf("%.4f ", df[:Wpwp_FRACT3][irow]))
        write(fid, @sprintf("%.4f ", df[:rough][irow]))
        write(fid, @sprintf("%.4f ", df[:snow_rough][irow]))
        write(fid, @sprintf("%.4f ", df[:annual_prec][irow]))
        write(fid, @sprintf("%.4f ", df[:resid_moist1][irow]))
        write(fid, @sprintf("%.4f ", df[:resid_moist2][irow]))
        write(fid, @sprintf("%.4f ", df[:resid_moist3][irow]))
        write(fid, @sprintf("%.4f \n", df[:fs_active][irow]))

    end

    close(fid)

    return nothing

end

"""
Read flux files (ascii outputs).
"""
function read_fluxes(path_sim, file)

    df = CSV.File(joinpath(path_sim, file), separator = '\t', header = false) |> DataFrame

    # Water balance mode on daily time step

    ncols = size(df, 2)

    if ncols == 25

        colnames = [:year, :month, :day, :prec, :evap,
        :runoff, :baseflow, :wdew, :soil_liq1, :soil_liq2,
        :soil_liq3, :net_short, :r_net, :evap_canop, :transp_veg,
        :evap_bare, :sub_canop, :sub_snow, :aero_resist, :surf_temp,
        :albedo, :rel_humid, :in_long, :air_temp, :wind]

        names!(df, colnames)

        df[:time] = DateTime.(df[:year], df[:month], df[:day])

    end

    # Water balance mode on sub-daily time step

    if ncols == 26

        colnames = [:year, :month, :day, :hour, :prec, :evap,
        :runoff, :baseflow, :wdew, :soil_liq1, :soil_liq2,
        :soil_liq3, :net_short, :r_net, :evap_canop, :transp_veg,
        :evap_bare, :sub_canop, :sub_snow, :aero_resist, :surf_temp,
        :albedo, :rel_humid, :in_long, :air_temp, :wind]

        names!(df, colnames)

        df[:time] = DateTime.(df[:year], df[:month], df[:day], df[:hour])

    end

    # Energy balance mode on daily time step

    if ncols == 31

        colnames = [:year, :month, :day, :prec, :evap,
        :runoff, :baseflow, :wdew, :soil_liq1, :soil_liq2,
        :soil_liq3, :rad_temp, :net_short, :r_net, :latent,
        :evap_canop, :transp_veg, :evap_bare, :sub_canop, :sub_snow,
        :sensible, :grnd_flux, :deltah, :fusion, :aero_resist,
        :surf_temp, :albedo, :rel_humid, :in_long, :air_temp,
        :wind]

        names!(df, colnames)

        df[:time] = DateTime.(df[:year], df[:month], df[:day])

    end

    # Energy balance mode on sub-daily time step

    if ncols == 32

        colnames = [:year, :month, :day, :hour, :prec, :evap,
        :runoff, :baseflow, :wdew, :soil_liq1, :soil_liq2,
        :soil_liq3, :rad_temp, :net_short, :r_net, :latent,
        :evap_canop, :transp_veg, :evap_bare, :sub_canop, :sub_snow,
        :sensible, :grnd_flux, :deltah, :fusion, :aero_resist,
        :surf_temp, :albedo, :rel_humid, :in_long, :air_temp,
        :wind]

        names!(df, colnames)

        df[:time] = DateTime.(df[:year], df[:month], df[:day], df[:hour])

    end

    return df

end

"""
Read snow files (ascii outputs).
"""
function read_snow(path_sim, file)

    df = CSV.File(joinpath(path_sim, file), separator = '\t', header = false) |> DataFrame

    # Water balance mode on daily time step

    ncols = size(df, 2)

    if ncols == 7

        colnames = [:year, :month, :day, :swe, :snow_depth, :snow_canopy, :snow_cover]

        names!(df, colnames)

        df[:time] = DateTime.(df[:year], df[:month], df[:day])

    end

    # Water balance mode on sub-daily time step

    if ncols == 8

        colnames = [:year, :month, :day, :hour, :swe, :snow_depth, :snow_canopy, :snow_cover]

        names!(df, colnames)

        df[:time] = DateTime.(df[:year], df[:month], df[:day], df[:hour])

    end

    # Energy balance mode on daily time step

    if ncols == 17

        colnames = [:year, :month, :day, :swe, :snow_depth, :snow_canopy, :snow_cover,
                    :advection, :deltacc, :snow_flux, :rfrz_energy, :melt_energy, :adv_sens, :latent_sub,
                    :snow_surf_temp, :snow_pack_temp, :snow_melt]

        names!(df, colnames)

        df[:time] = DateTime.(df[:year], df[:month], df[:day])

    end

    # Energy balance mode on sub-daily time step

    if ncols == 18

        colnames = [:year, :month, :day, :hour, :swe, :snow_depth, :snow_canopy, :snow_cover,
                    :advection, :deltacc, :snow_flux, :rfrz_energy, :melt_energy, :adv_sens, :latent_sub,
                    :snow_surf_temp, :snow_pack_temp, :snow_melt]

        names!(df, colnames)

        df[:time] = DateTime.(df[:year], df[:month], df[:day], df[:hour])

    end

    return df

end

"""
Read vic results (binary outputs).
"""
function read_results(path, file, ntime, ndata)

    fid = open(joinpath(path, file), "r")

    time_array = Int32[]
    data_array = Float64[]

    while !eof(fid)

        for i in 1:ntime
            push!(time_array, read(fid, Int32))
        end

        for i in 1:ndata
            push!(data_array, read(fid, Float64))
        end

    end

    close(fid)

    return time_array, data_array

end

"""
Read vic forcings.
"""
function read_vic_forcing(file)

    prec = Float64[]
    tmin = Float64[]
    tmax = Float64[]
    wind = Float64[]

    fid = open(file, "r")

    while !eof(fid)

        push!(prec, convert(Float64, read(fid, Int16))/100.0)
        push!(tmin, convert(Float64, read(fid, Int16))/100.0)
        push!(tmax, convert(Float64, read(fid, Int16))/100.0)
        push!(wind, convert(Float64, read(fid, Int16))/100.0)

    end

    close(fid)

    return prec, tmin, tmax, wind

end

"""
Write vic forcings.
"""
function write_vic_forcing(file, prec, tmin, tmax, wind)

    fid = open(file, "w")

    for i in 1:length(prec)

        write(fid, round(Int16, 100*prec[i]))
        write(fid, round(Int16, 100*tmin[i]))
        write(fid, round(Int16, 100*tmax[i]))
        write(fid, round(Int16, 100*wind[i]))

    end

    close(fid)

    return nothing

end

"""
Read mtclim output.
"""
function read_mtclim_hbv(file)

    #=
    OUTVAR		OUT_PREC	    * OUT_TYPE_USINT	40      incoming precipitation [mm]
    OUTVAR		OUT_AIR_TEMP        * OUT_TYPE_SINT	100     air temperature [C]
    OUTVAR		OUT_SHORTWAVE	    * OUT_TYPE_USINT	50      incoming shortwave [W/m2]
    OUTVAR		OUT_LONGWAVE        * OUT_TYPE_USINT	80      incoming longwave [W/m2]
    OUTVAR		OUT_PRESSURE	    * OUT_TYPE_USINT	100     near surface atmospheric pressure [kPa]
    OUTVAR		OUT_QAIR	    * OUT_TYPE_USINT	100000  specific humidity [kg/kg]
    OUTVAR		OUT_VP		    * OUT_TYPE_USINT	100     near surface vapor pressure [kPa]
    OUTVAR		OUT_WIND	    * OUT_TYPE_USINT	100     near surface wind speed [m/s]
    =#

    prec = []
    tair = []
    iswr = []
    ilwr = []
    pres = []
    qair = []
    vp = []
    wind = []

    fid = open(file, "r")

    while !eof(fid)

        push!(prec , convert(Float64, read(fid, UInt16))/40.0)
        push!(tair , convert(Float64, read(fid, Int16))/100.0)
        push!(iswr , convert(Float64, read(fid, UInt16))/50.0)
        push!(ilwr , convert(Float64, read(fid, UInt16))/80.0)
        push!(pres , convert(Float64, read(fid, UInt16))/100.0)
        push!(qair , convert(Float64, read(fid, UInt16))/100000.0)
        push!(vp ,   convert(Float64, read(fid, UInt16))/100.0)
        push!(wind , convert(Float64, read(fid, UInt16))/100.0)

    end

    close(fid)

    return prec, tair, iswr, ilwr, pres, qair, vp, wind

end

"""
Read mtclim output.
"""
function read_mtclim_fsm(file)

    #=
    OUTVAR		OUT_RAINF	* OUT_TYPE_USINT	40    [mm]
    OUTVAR		OUT_SNOWF	* OUT_TYPE_USINT	40    [mm]
    OUTVAR		OUT_AIR_TEMP	* OUT_TYPE_SINT		100   [C]
    OUTVAR		OUT_SHORTWAVE	* OUT_TYPE_USINT	50    [W/m2]
    OUTVAR		OUT_LONGWAVE	* OUT_TYPE_USINT	80    [W/m2]
    OUTVAR		OUT_PRESSURE	* OUT_TYPE_USINT	100   [kPa]
    OUTVAR		OUT_REL_HUMID	* OUT_TYPE_USINT	100   [%]
    OUTVAR		OUT_WIND	* OUT_TYPE_USINT	100   [m/s]
    =#

    rainf = Float64[]
    snowf = Float64[]
    tair = Float64[]
    iswr = Float64[]
    ilwr = Float64[]
    pres = Float64[]
    rhum = Float64[]
    wind = Float64[]

    fid = open(file, "r")

    while !eof(fid)

        push!(rainf , convert(Float64, read(fid, UInt16))/40.0)
        push!(snowf , convert(Float64, read(fid, UInt16))/40.0)
        push!(tair  , convert(Float64, read(fid, Int16))/100.0)
        push!(iswr  , convert(Float64, read(fid, UInt16))/50.0)
        push!(ilwr  , convert(Float64, read(fid, UInt16))/80.0)
        push!(pres  , convert(Float64, read(fid, UInt16))/100.0)
        push!(rhum  , convert(Float64, read(fid, UInt16))/100.0)
        push!(wind  , convert(Float64, read(fid, UInt16))/100.0)

    end

    close(fid)

    return rainf, snowf, tair, iswr, ilwr, pres, rhum, wind

end

"""
Read vegetation parameter file and return as dataframe.
"""
function read_veg_param(path)

    # Standard file name

    file = "params/veg_param"

    # Empty dataframe

    df = DataFrame(gridcel = [],
                   Nveg = [],
                   veg_class = [],
                   Cv = [],
                   root_depth1 = [],
                   root_depth2 = [],
                   root_frac1 = [],
                   root_frac2 = [],
                   lai1 = [],
                   lai2 = [],
                   lai3 = [],
                   lai4 = [],
                   lai5 = [],
                   lai6 = [],
                   lai7 = [],
                   lai8 = [],
                   lai9 = [],
                   lai10 = [],
                   lai11 = [],
                   lai12 = [])

    # Read file

    str = readlines(joinpath(path, file))

    i = 1

    while i < length(str)

        # Read gridcel and number of vegetation classes

        tmp = parse.(split(str[i]))

        gridcel = tmp[1]
        Nveg = tmp[2]

        for j = 1:2:(2*Nveg-1)

            # Process row with vegetation class ...

            tmp = parse.(split(str[i+j]))

            veg_class = tmp[1]
            Cv = tmp[2]
            root_depth = tmp[3:4]
            root_frac = tmp[5:6]

            # Process row with leaf area index

            tmp = parse.(split(str[i+j+1]))

            lai = tmp

            # Add all together

            row = vcat(gridcel, Nveg, veg_class, Cv, root_depth, root_frac, lai)

            push!(df, row)

        end

        i += 2*Nveg + 1

    end

    return df

end

"""
Write vegetation parameter file given as a dataframe.
"""
function write_veg_params(path, df)

    fid = open(joinpath(path, "params/veg_param"), "w")

    for gridcel in unique(df[:gridcel])

        irow = findin(df[:gridcel], gridcel)

        df_tmp = df[irow, :]

        Nveg = df_tmp[:Nveg][1]

        write(fid, @sprintf("%.0f %.0f\n", gridcel, Nveg))

        for row in eachrow(df_tmp)

            write(fid, @sprintf("%6.0f %0.3f %0.2f %0.2f %0.2f %0.2f\n",
                                row[:veg_class], row[:Cv], row[:root_depth1], row[:root_depth2], row[:root_frac1], row[:root_frac2]))

            write(fid, @sprintf("       %-0.3f %0.3f %0.3f %0.3f %0.3f %0.3f %0.3f %0.3f %0.3f %0.3f %0.3f %0.3f\n",
                                row[:lai1], row[:lai2], row[:lai3], row[:lai4], row[:lai5], row[:lai6],
                                row[:lai7], row[:lai8], row[:lai9], row[:lai10], row[:lai11], row[:lai12]))

        end

    end

    close(fid)

end

"""
Read all snow files in a directory (ascii files).
"""
function read_all_snow(path)

    files = readdir(path)

    files = filter(x->contains(x, "snow"), files)

    df_tmp = read_snow(path, files[1])

    ifirst = find(names(df_tmp) .== :swe)[1]
    ilast = find(names(df_tmp) .== :snow_melt)[1]

    data = convert(Array, df_tmp[ifirst:ilast])/length(files)

    for i in 2:length(files)
        @show i
        df_tmp = read_snow(path, files[i])
        data += convert(Array, df_tmp[ifirst:ilast])/length(files)
    end

    df = DataFrame(data)

    names!(df, names(df_tmp[ifirst:ilast]))

    df[:time] = df_tmp[:time]

    return df

end

"""
Read all flux files in a directory (ascii files).
"""
function read_all_fluxes(path)

    files = readdir(path)

    files = filter(x->contains(x, "fluxes"), files)

    df_tmp = read_fluxes(path, files[1])

    ifirst = find(names(df_tmp) .== :prec)[1]
    ilast = find(names(df_tmp) .== :wind)[1]

    data = convert(Array, df_tmp[ifirst:ilast])/length(files)

    for i in 2:length(files)
        @show i
        df_tmp = read_fluxes(path, files[i])
        data += convert(Array, df_tmp[ifirst:ilast])/length(files)
    end

    df = DataFrame(data)

    names!(df, names(df_tmp[ifirst:ilast]))

    df[:time] = df_tmp[:time]

    return df

end

#=
"""
Read all result files in a directory (binary files).
"""
function read_all_results(path)

    info("Processing data in folder: $(path)")

    # Data info

    colnames = [:prec, :evap, :runoff, :baseflow, :wdew, :net_short, :r_net,
                :latent, :evap_canop, :transp_veg, :evap_bare, :sub_canop,
                :sub_snow, :rel_humid, :in_long, :wind, :swe, :snow_depth,
                :snow_canopy, :snow_cover]

    # Number of time and data columns in results files

    ntime = 3

    ndata = length(colnames)

    # Read files

    files = readdir(path)

    time, data = read_results(path, files[1], ntime, ndata)

    data = data / length(files)

    @showprogress 1 "Processing..." for i in 2:length(files)
        time_tmp, data_tmp = read_results(path, files[i], ntime, ndata)
        data += data_tmp / length(files)
    end

    # Reshape vectors

    time = reshape(time, ntime, :)'
    data = reshape(data, ndata, :)'

    # Create data frame

    df = DataFrame(data)

    names!(df, colnames)

    df[:time] = DateTime.(time[:,1], time[:,2], time[:,3])

    return df

end
=#

"""
Read all result files in a directory (binary files).
"""
function read_all_results(path, watershed, resolution, wsh_name)

    info("Processing data in folder: $(path)")

    # Data info

    colnames = [:prec, :air_temp, :rainf, :snowf, :wind, :runoff, :baseflow, :inflow,
                :evap, :sub_snow, :transp_veg, :sub_canop, :zwt_lumped, :soil_wet, :swe,
                :snow_cover, :snow_canopy, :wdew, :latent, :sensible, :net_long, :net_short,
                :pet_h2osurf, :pet_short, :pet_natveg, :soil_liq1, :soil_liq2, :soil_liq3,
                :evap_canop, :evap_bare]

    # Number of time and data columns in results files

    ntime = 3

    ndata = length(colnames)

    # Read soil parameter file

    soil_param = read_soil_params(joinpath(path, watershed, resolution))

    lat = soil_param[:lat]
    lon = soil_param[:lon]
    elev = soil_param[:elev]
    gridcel = soil_param[:gridcel]
    name = wsh_name[watershed]

    # Load one file and allocate output array

    lat_str = @sprintf("%0.5f", lat[1])
    lon_str = @sprintf("%0.5f", lon[1])

    file_res = joinpath(path, watershed, resolution, "results/results_$(lat_str)_$(lon_str)")

    time_raw, data_raw = read_results(path, file_res, ntime, ndata)

    data = reshape(data_raw, ndata, :)'

    data_all = zeros(size(data,1), size(data,2), size(soil_param,1))

    # Load results from all files

    @showprogress 1 "Processing..." for i in 1:size(soil_param,1)

        lat_str = @sprintf("%0.5f", lat[i])
        lon_str = @sprintf("%0.5f", lon[i])

        file_res = joinpath(path, watershed, resolution, "results/results_$(lat_str)_$(lon_str)")

        time_raw, data_raw = read_results(path, file_res, ntime, ndata)

        data = reshape(data_raw, ndata, :)'

        data_all[:,:,i] = data

    end

    # Average data spatially

    data_mean = mean(data_all, 3)[:,:,1]

    # Handle time

    time = reshape(time_raw, ntime, :)'

    time = DateTime.(time[:,1], time[:,2], time[:,3])

    # Create data frame

    df_vic = DataFrame(data_mean)

    names!(df_vic, colnames)

    df_vic[:time] = time

    # Create result struct

    res_vic = VicRes(name, watershed, lat, lon, gridcel, elev, data_all, data_mean, colnames, time)

    return df_vic, res_vic

end


text:
using Flux.Optimise: update!
using Statistics: std, mean

function policyupdate!(rewards, log_probs)
    discounted_rewards = Float32[calcreward(rewards[t:end]) for t=1:length(rewards)]

    # Using normalized discounted rewards, because Andrej Karpathy said so.
    # This is a hack to control the variance of the policy gradient estimator.
    norm_discounted_rewards = (discounted_rewards .- mean(discounted_rewards)) ./ (std(discounted_rewards) .+ 1f-6)

    policy_grads = -log_probs .* norm_discounted_rewards

    ps = params(model)
    gs = gradient(()->sum(policy_grads), ps)

    update!(opt, ps, gs)
end

function calcreward(rewards_subset)
    γ_subarray = ones(Float32, length(rewards_subset)-1)
    γ_subarray *= γ
    γ_array = cumprod(vcat(Float32[1], γ_subarray))
    discouonted_goal = sum(γ_array .* rewards_subset)
    return discouonted_goal
end

ext: jl


text:
using ArviZPlots
using Test

@testset "ArviZPlots.jl" begin
    # Write your tests here.
end


text:
module Interpolation

    using ..Utils

    export Interpolator, HermiteInterpolation

    include("interpolation/interpolation.jl")
    include("interpolation/hermite_interpolation.jl")

end


text:
# ### m8.1t

using Pkg, DrWatson

@quickactivate "StatisticalRethinkingTuring"
using Turing
using StatisticalRethinking
Turing.setprogress!(false)

@model ppl8_2(y) = begin

    σ ~ FlatPos(0.0) # improper prior with probability one everywhere above 0.0
    α ~ Flat() # improper prior with pobability one everywhere

    y .~ Normal(α, σ)
end

y = [-1,1];

# Sample

m8_2t = ppl8_2(y)
nchains = 4; sampler = NUTS(0.65); nsamples=2000
chns8_2t = mapreduce(c -> sample(m8_2t, sampler, nsamples), chainscat, 1:nchains)

# End of `08/m8.2t.jl`


text:
"""
    Vocabulary(word,unk_cutoff =1 ,unk_label = "<unk>")

Stores language model vocabulary.
Satisfies two common language modeling requirements for a vocabulary:
- When checking membership and calculating its size, filters items
by comparing their counts to a cutoff value.
Adds a special "unknown" token which unseen words are mapped to.

# Example
```julia-repl
julia> words = ["a", "c", "-", "d", "c", "a", "b", "r", "a", "c", "d"]
julia> vocabulary = Vocabulary(words, 2)
  Vocabulary(Dict("<unk>"=>1,"c"=>3,"a"=>3,"d"=>2), 2, "<unk>")

julia> vocabulary.vocab
  Dict{String,Int64} with 4 entries:
   "<unk>" => 1
   "c"     => 3
   "a"     => 3
   "d"     => 2

Tokens with counts greater than or equal to the cutoff value will
be considered part of the vocabulary.
julia> vocabulary.vocab["c"]
 3

julia> "c" in keys(vocabulary.vocab)
 true

julia> vocabulary.vocab["d"]
 2

julia> "d" in keys(vocabulary.vocab)
 true

Tokens with frequency counts less than the cutoff value will be considered not
part of the vocabulary even though their entries in the count dictionary are
preserved.
julia> "b" in keys(vocabulary.vocab)
 false

julia> "<unk>" in keys(vocabulary.vocab)
 true

We can look up words in a vocabulary using its `lookup` method.
"Unseen" words (with counts less than cutoff) are looked up as the unknown label.
If given one word (a string) as an input, this method will return a string.
julia> lookup("a")
 'a'

julia> word = ["a", "-", "d", "c", "a"]

julia> lookup(vocabulary ,word)
 5-element Array{Any,1}:
  "a"
  "<unk>"
  "d"
  "c"
  "a"

If given a sequence, it will return an Array{Any,1} of the looked up words as shown above.

It's possible to update the counts after the vocabulary has been created.
julia> update(vocabulary,["b","c","c"])
 1

julia> vocabulary.vocab["b"]
 1
```
"""
mutable struct Vocabulary
vocab::Dict{String, Int64}
unk_cutoff::Int
unk_label::String
allword::Array{String, 1}
end

function Vocabulary(word::Vector{T}, unk_cutoff=1, unk_label="<unk>") where { T <: AbstractString}
    if unk_label in word
        error("unk_label is in vocab")
    else
    word= push!(word, unk_label)
    end
    vocab = countmap(word)
    for value in vocab
        if value[2]<unk_cutoff && value[1] != unk_label
            delete!(vocab, value[1])
        end
    end
    Vocabulary(vocab, unk_cutoff, unk_label, word)
end

function update(vocab::Vocabulary, words)
    vocab.allword = append!(vocab.allword, words)
    vocab.vocab=addcounts!(vocab.vocab, words)
end

"""
lookup a sequence or words in the vocabulary

Return an Array of String
"""
function lookup(voc::Vocabulary, word::Vector{T}) where { T <: AbstractString}
    look = []
    for w in word
        if w in keys(voc.vocab)
            push!(look, w)
        else
            push!(look, voc.unk_label)
        end
    end
    return look
end


text:
show(io::IO, c::Colorant)              = _show(io, c)
show(io::IO, c::ColorantNormed)        = show_normed(io, c)
showcompact(io::IO, c::Colorant)       = _showcompact(io, c)
showcompact(io::IO, c::ColorantNormed) = show_normed(io, c)

for N = 1:4
    component = N >= 3 ? (:comp1, :comp2, :comp3, :alpha) : (:comp1, :alpha)
    printargs = Array{Any}(2, N)
    for i = 1:N
        printargs[1,i] = :(show(io, $(component[i])(c)))
        chr = i < N ? ',' : ')'
        printargs[2,i] = :(print(io, $chr))
    end
    @eval begin
        function _show{T}(io::IO, c::Colorant{T,$N})
            print(io, colorant_string(typeof(c)), "{", T, "}(")
            $(printargs[:]...)
        end
    end
    for i = 1:N
        printargs[1,i] = :(showcompact(io, $(component[i])(c)))
    end
    @eval begin
        function _showcompact{T}(io::IO, c::Colorant{T,$N})
            print(io, colorant_string(typeof(c)), "{", T, "}(")
            $(printargs[:]...)
        end
        # Special handling for Normed types: don't print the giant type name
        function show_normed{T,f}(io::IO, c::Colorant{FixedPointNumbers.Normed{T,f},$N})
            print(io, colorant_string(typeof(c)), '{')
            FixedPointNumbers.showtype(io, eltype(typeof(c)))
            print(io, "}(")
            $(printargs[:]...)
        end
    end
end


text:
"""
    Model Data
"""

struct DynamicsDerivativesData{X, U, W}
    fx::Vector{X}
    fu::Vector{U}
	fw::Vector{W}
end

function dynamics_derivatives_data(model::Model, T;
	n = [model.n for t = 1:T],
	m = [model.m for t = 1:T-1],
	d = [model.d for t = 1:T-1])

	fx = [zeros(n[t + 1], n[t]) for t = 1:T-1]
    fu = [zeros(n[t + 1], m[t]) for t = 1:T-1]
	fw = [zeros(n[t + 1], d[t]) for t = 1:T-1]

    DynamicsDerivativesData(fx, fu, fw)
end

struct ObjectiveDerivativesData{X, U, XX, UU, UX}
    gx::Vector{X}
    gu::Vector{U}
    gxx::Vector{XX}
    guu::Vector{UU}
    gux::Vector{UX}
end

function objective_derivatives_data(model::Model, T;
	n = [model.n for t = 1:T],
	m = [model.m for t = 1:T-1])

	gx = [ones(n[t]) for t = 1:T]
    gu = [ones(m[t]) for t = 1:T-1]
    gxx = [ones(n[t], n[t]) for t = 1:T]
    guu = [ones(m[t], m[t]) for t = 1:T-1]
    gux = [ones(m[t], n[t]) for t = 1:T-1]

    ObjectiveDerivativesData(gx, gu, gxx, guu, gux)
end

struct ModelData{X, U, D, S}
    # current trajectory
    x::Vector{X}
    u::Vector{U}

    # disturbance trajectory
    w::Vector{D}

    # time step
    h::S

    # horizon
    T::Int

    # nominal trajectory
    x̄::Vector{X}
    ū::Vector{U}

    # dynamics model
    model::Model
	n::Vector{Int}
	m::Vector{Int}
	d::Vector{Int}

    # objective
    obj::Objective

    # dynamics derivatives data
    dyn_deriv::DynamicsDerivativesData

    # objective derivatives data
    obj_deriv::ObjectiveDerivativesData

    # z = (x1...,xT,u1,...,uT-1) | Δz = (Δx1...,ΔxT,Δu1,...,ΔuT-1)
    z::Vector{S}

	analytical_derivatives::Bool
end

ModelsData = Vector{ModelData}

function model_data(model, obj, w, h, T;
	n = [model.n for t = 1:T],
	m = [model.m for t = 1:T-1],
	d = [model.d for t = 1:T-1],
	analytical_derivatives = false)

    num_var = sum(n) + sum(m)

	x = [zeros(n[t]) for t = 1:T]
    u = [zeros(m[t]) for t = 1:T-1]

    x̄ = [zeros(n[t]) for t = 1:T]
    ū = [zeros(m[t]) for t = 1:T-1]

    dyn_deriv = dynamics_derivatives_data(model, T, n = n, m = m, d = d)
    obj_deriv = objective_derivatives_data(model, T, n = n, m = m)

    z = zeros(num_var)

    ModelData(x, u, w, h, T, x̄, ū, model, n, m, d, obj, dyn_deriv, obj_deriv, z,
		analytical_derivatives)
end

# function Δz!(m_data::ModelData)
# 	n = m_data.n
# 	m = m_data.m
# 	T = m_data.T
#
#     for t = 1:T
#         idx_x = (t == 1 ? 0 : (t - 1) * n[t-1]) .+ (1:n[t])
#         m_data.z[idx_x] = m_data.x[t] - m_data.x̄[t]
#
#         t == T && continue
#
#         idx_u = sum(n) + (t == 1 ? 0 : (t - 1) * m[t-1]) .+ (1:m[t])
#         m_data.z[idx_u] = m_data.u[t] - m_data.ū[t]
#     end
# end

function objective(data::ModelData; mode = :nominal)
    if mode == :nominal
        return objective(data.obj, data.x̄, data.ū)
    elseif mode == :current
        return objective(data.obj, data.x, data.u)
    end
end

"""
    Policy Data
"""
struct PolicyData{N, M, NN, MM, MN, NNN, MNN}
    # policy
    K::Vector{MN}
    k::Vector{M}

    K_cand::Vector{MN}
    k_cand::Vector{M}

    # value function approximation
    P::Vector{NN}
    p::Vector{N}

    # state-action value function approximation
    Qx::Vector{N}
    Qu::Vector{M}
    Qxx::Vector{NN}
    Quu::Vector{MM}
    Qux::Vector{MN}

	xx̂_tmp::Vector{NNN}
	ux̂_tmp::Vector{MNN}
	uu_tmp::Vector{MM}
	ux_tmp::Vector{MN}
end

function policy_data(model::Model, T;
	n = [model.n for t = 1:T],
	m = [model.m for t = 1:T-1])

	K = [zeros(m[t], n[t]) for t = 1:T-1]
    k = [zeros(m[t]) for t = 1:T-1]

    K_cand = [zeros(m[t], n[t]) for t = 1:T-1]
    k_cand = [zeros(m[t]) for t = 1:T-1]

    P = [zeros(n[t], n[t]) for t = 1:T]
    p = [zeros(n[t]) for t = 1:T]

    Qx = [zeros(n[t]) for t = 1:T-1]
    Qu = [zeros(m[t]) for t = 1:T-1]
    Qxx = [zeros(n[t], n[t]) for t = 1:T-1]
    Quu = [zeros(m[t], m[t]) for t = 1:T-1]
    Qux = [zeros(m[t], n[t]) for t = 1:T-1]

	xx̂_tmp = [zeros(n[t], n[t + 1]) for t = 1:T-1]
	ux̂_tmp = [zeros(m[t], n[t + 1]) for t = 1:T-1]
	uu_tmp = [zeros(m[t], m[t]) for t = 1:T-1]
	ux_tmp = [zeros(m[t], n[t]) for t = 1:T-1]

    PolicyData(K, k, K_cand, k_cand, P, p, Qx, Qu, Qxx, Quu, Qux, xx̂_tmp, ux̂_tmp, uu_tmp, ux_tmp)
end

"""
    Solver Data
"""
mutable struct SolverData{T}
    obj::T              # objective value
    gradient::Vector{T} # Lagrangian gradient
	c_max::T            # maximum constraint violation

    idx_x::Vector       # indices for state trajectory
    idx_u::Vector       # indices for control trajectory

    α::T                # step length
    status::Bool        # solver status

	cache::Dict         #
end

cache = Dict(:obj => [], :gradient => [], :c_max => [], :α => [])

function solver_data(model::Model, T;
	n = [model.n for t = 1:T],
	m = [model.m for t = 1:T-1])

    num_var = sum(n) + sum(m)

    idx_x = [sum(n[1:(t-1)]) .+ (1:n[t]) for t = 1:T]
    idx_u = [sum(n) + sum(m[1:(t-1)]) .+ (1:m[t]) for t = 1:T-1]

    obj = Inf
	c_max = 0.0
    gradient = zeros(num_var)
	cache = Dict(:obj => [], :gradient => [], :c_max => [], :α => [])

    SolverData(obj, gradient, c_max, idx_x, idx_u, 1.0, false, cache)
end

function cache!(data::SolverData)
	push!(data.cache[:obj], data.obj)
	push!(data.cache[:gradient], data.gradient)
	push!(data.cache[:c_max], data.c_max)
	push!(data.cache[:α], data.α)
end

function objective!(s_data::SolverData, m_data::ModelData; mode = :nominal)
	if mode == :nominal
		s_data.obj = objective(m_data.obj, m_data.x̄, m_data.ū)
	elseif mode == :current
		s_data.obj = objective(m_data.obj, m_data.x, m_data.u)
	end

	if m_data.obj isa AugmentedLagrangianCosts
		s_data.c_max = constraint_violation(m_data.obj.cons,
			m_data.x, m_data.u,
			norm_type = Inf)
	end

	return s_data.obj
end

function Δz!(m_data::ModelData, p_data::PolicyData, s_data::SolverData)
	n = m_data.n
	m = m_data.m
	T = m_data.T

    m_data.z[s_data.idx_x[1]] .= 0.0

    for t = 1:T-1
        m_data.z[s_data.idx_u[t]] .= p_data.K[t] * m_data.z[s_data.idx_x[t]] + p_data.k[t]
        m_data.z[s_data.idx_x[t+1]] .= m_data.dyn_deriv.fx[t] * m_data.z[s_data.idx_x[t]] + m_data.dyn_deriv.fu[t] * m_data.z[s_data.idx_u[t]]
    end
end

"""
    Problem Data
"""
struct ProblemData
	p_data
	m_data
	s_data
end

function problem_data(model::Model, obj::StageCosts, x̄, ū, w, h, T;
	n = [model.n for t = 1:T],
	m = [model.m for t = 1:T-1],
	d = [model.d for t = 1:T-1],
	analytical_dynamics_derivatives = false)

	# allocate policy data
    p_data = policy_data(model, T, n = n, m = m)

    # allocate model data
    m_data = model_data(model, obj, w, h, T, n = n, m = m, d = d,
		analytical_derivatives = analytical_dynamics_derivatives)
    m_data.x̄ .= x̄
    m_data.ū .= ū

    # allocate solver data
    s_data = solver_data(model, T, n = n, m = m)

	ProblemData(p_data, m_data, s_data)
end

function problem_data(m_data::ModelsData;
	n = m_data[1].n,
	m = m_data[1].m)

	model = m_data[1].model
	T = m_data[1].T

	# allocate policy data
    p_data = policy_data(model, T, n = n, m = m)

    # allocate solver data
    s_data = solver_data(model, T, n = n, m = m)

	ProblemData(p_data, m_data, s_data)
end

function nominal_trajectory(prob::ProblemData)
	return prob.m_data.x̄, prob.m_data.ū
end

function current_trajectory(prob::ProblemData)
	return prob.m_data.x, prob.m_data.u
end

function problem_data(model, obj::StageCosts, con_set::ConstraintSet,
		x̄, ū, w, h, T;
		n = [model.n for t = 1:T],
		m = [model.m for t = 1:T-1],
		d = [model.d for t = 1:T-1],
		analytical_dynamics_derivatives = false)

	# constraints
	c_data = constraints_data(model, [c.p for c in con_set], T, n = n, m = m)
	cons = StageConstraints(con_set, c_data, T)

	# augmented Lagrangian
	obj_al = augmented_lagrangian(obj, cons)

	# allocate policy data
    p_data = policy_data(model, T, n = n, m = m)

    # allocate model data
    m_data = model_data(model, obj_al, w, h, T, n = n, m = m, d = d,
		analytical_derivatives = analytical_dynamics_derivatives)
    m_data.x̄ .= x̄
    m_data.ū .= ū

    # allocate solver data
    s_data = solver_data(model, T, n = n, m = m)

	ProblemData(p_data, m_data, s_data)
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
function capital_grid_test()
    kGrid = init_test_capital_grid();
    @test validate_capital_grid(kGrid);

    t = 5;
    kGridV = make_k_grid(kGrid, t);
    @test validate_capital_grid(kGridV);
end

@testset "Capital grid" begin
    capital_grid_test();
end

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


text:
using Luxor, FileIO, ImageCore, Statistics, DelimitedFiles, Serialization
fontheight = 80

for font in ["Courier"]
    @info "Rendering $font characters"
    fontdir = joinpath(@__DIR__,"fonts",font)
    !isdir(fontdir) && mkdir(fontdir)
    xFirsts = Int[]
    xLasts = Int[]
    yFirsts = Int[]
    yLasts = Int[]
    for charIDX in 0:255
        char = Char(charIDX)
        file = joinpath(fontdir,"$charIDX.png")
        Drawing(150, 150, file)
        origin()
        background("black")
        fontface(font)
        setcolor("white")
        setopacity(1)
        fontsize(fontheight)
        text(string(char), Point(0,0), halign = :bottom,  valign = :left)
        finish()
        img = load(file)
        pixcol = mean(img, dims=1)[1,:]
        xOccupied = findall(pixcol .!= RGB(0,0,0))
        if length(xOccupied)>2
            push!(xFirsts, Int(xOccupied[1]))
            push!(xLasts, Int(xOccupied[end]))
        end
        pixrow = mean(img, dims=2)[:,1]
        yOccupied = findall(pixrow .!= RGB(0,0,0))
        if length(yOccupied)>2
            push!(yFirsts, Int(yOccupied[1]))
            push!(yLasts, Int(yOccupied[end]))
        end
    end
    x_first_pix = minimum(xFirsts)
    x_last_pix = maximum(xLasts)
    y_first_pix = minimum(yFirsts)
    y_last_pix = maximum(yLasts)

    charIDX = 0
    charStack = Vector{Array{Float64}}(undef,256)

    for charIDX in 0:255
        file = joinpath(fontdir,"$charIDX.png")
        img = load(file)
        charStack[charIDX+1] = round.(UInt8,channelview(img[y_first_pix:y_last_pix, x_first_pix:x_last_pix])[1,:,:] .* 255)
        rm(file)
    end
    Serialization.serialize(joinpath(fontdir,"charStack"), charStack)
end


text:
module WignerSymbols
export δ, Δ, clebschgordan, wigner3j, wigner6j, racahV, racahW, HalfInteger

using HalfIntegers
using RationalRoots
using LRUCache
const RRBig = RationalRoot{BigInt}
import RationalRoots: _convert

include("growinglist.jl")
include("bigint.jl") # additional GMP BigInt functionality not wrapped in Base.GMP.MPZ
include("primefactorization.jl")
convert(BigInt, primefactorial(401)) # trigger compilation and generate some fixed data

const Key3j = Tuple{UInt,UInt,UInt,Int,Int}
const Key6j = NTuple{6,UInt}

const Wigner3j = LRU{Key3j,Tuple{Rational{BigInt},Rational{BigInt}}}(; maxsize = 10^6)
const Wigner6j = LRU{Key6j,Tuple{Rational{BigInt},Rational{BigInt}}}(; maxsize = 10^6)

function set_buffer3j_size!(; maxsize)
    resize!(Wigner3j; maxsize = maxsize)
end
function set_buffer6j_size!(; maxsize)
    resize!(Wigner6j; maxsize = maxsize)
end

# check integerness and correctness of (j,m) angular momentum
ϵ(j, m) = (abs(m) <= j && ishalfinteger(j) && isinteger(j-m) && isinteger(j+m))

# check triangle condition
"""
    δ(j₁, j₂, j₃) -> ::Bool

Checks the triangle conditions `j₃ <= j₁ + j₂`, `j₁ <= j₂ + j₃` and `j₂ <= j₃ + j₁`.
"""
δ(j₁, j₂, j₃) = (j₃ <= j₁ + j₂) && (j₁ <= j₂ + j₃) && (j₂ <= j₃ + j₁) && isinteger(j₁+j₂+j₃)

# triangle coefficient
"""
    Δ(T::Type{<:Real} = RationalRoot{BigInt}, j₁, j₂, j₃) -> ::T

Computes the triangle coefficient
`Δ(j₁, j₂, j₃) = √((j₁+j₂-j₃)!*(j₁-j₂+j₃)!*(j₂+j₃-j₁)! / (j₁+j₂+j₃+1)!)`
as a type `T` real number.

Returns `zero(T)` if the triangle condition `δ(j₁, j₂, j₃)` is not satisfied, but
throws a `DomainError` if the `jᵢ`s are not (half)integer
"""
Δ(j₁, j₂, j₃) = Δ(RRBig, j₁, j₂, j₃)
function Δ(T::Type{<:Real}, j₁, j₂, j₃)
    for jᵢ in (j₁, j₂, j₃)
        (ishalfinteger(jᵢ) && jᵢ >= zero(jᵢ)) || throw(DomainError("invalid jᵢ", jᵢ))
    end
    if !δ(j₁, j₂, j₃)
        return zero(T)
    end
    n, d = Δ²(j₁, j₂, j₃)
    r = Base.unsafe_rational(n, d)
    return convert(T, signedroot(RationalRoot{BigInt}, r))
end

"""
    wigner3j(T::Type{<:Real} = RationalRoot{BigInt}, j₁, j₂, j₃, m₁, m₂, m₃ = -m₂-m₁) -> ::T

Compute the value of the Wigner-3j symbol
    ⎛ j₁  j₂  j₃ ⎞
    ⎝ m₁  m₂  m₃ ⎠
as a type `T` real number. By default, `T = RationalRoot{BigInt}` and `m₃ = -m₁-m₂`.

Returns `zero(T)` if the triangle condition `δ(j₁, j₂, j₃)` is not satisfied, but
throws a `DomainError` if the `jᵢ`s and `mᵢ`s are not (half)integer or `abs(mᵢ) > jᵢ`.
"""
wigner3j(j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m₂) = wigner3j(RRBig, j₁, j₂, j₃, m₁, m₂, m₃)
function wigner3j(T::Type{<:Real}, j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m₂)
    # check angular momenta
    for (jᵢ,mᵢ) in ((j₁, m₁), (j₂, m₂), (j₃, m₃))
        ϵ(jᵢ, mᵢ) || throw(DomainError((jᵢ, mᵢ), "invalid combination (jᵢ, mᵢ)"))
    end
    return _wigner3j(T, HalfInteger.((j₁, j₂, j₃, m₁, m₂, m₃))...)
end

function _wigner3j(T::Type{<:Real}, j₁::HalfInteger, j₂::HalfInteger, j₃::HalfInteger,
                                    m₁::HalfInteger, m₂::HalfInteger, m₃::HalfInteger)
    # check triangle condition and m₁+m₂+m₃ == 0
    if !δ(j₁, j₂, j₃) || !iszero(m₁+m₂+m₃)
        return zero(T)
    end

    # we reorder such that j₁ >= j₂ >= j₃ and m₁ >= 0 or m₁ == 0 && m₂ >= 0
    j₁, j₂, j₃, m₁, m₂, m₃, sgn = reorder3j(j₁, j₂, j₃, m₁, m₂, m₃)
    # TODO: do we also want to use Regge symmetries?
    α₁ = convert(Int, j₂ - m₁ - j₃ ) # can be negative
    α₂ = convert(Int, j₁ + m₂ - j₃ ) # can be negative
    β₁ = convert(UInt, j₁ + j₂ - j₃ )
    β₂ = convert(UInt, j₁ - m₁ )
    β₃ = convert(UInt, j₂ + m₂ )

    # extra sign in definition: α₁ - α₂ = j₁ + m₂ - j₂ + m₁ = j₁ - j₂ + m₃
    sgn = isodd(α₁ - α₂) ? -sgn : sgn

    # dictionary lookup or compute
    if haskey(Wigner3j, (β₁, β₂, β₃, α₁, α₂))
        r, s = Wigner3j[(β₁, β₂, β₃, α₁, α₂)]
    else
        s1n, s1d = Δ²(j₁, j₂, j₃)
        s2n = prod(map(primefactorial, (β₂, β₁ - α₁, β₁ - α₂, β₃, β₃ - α₁, β₂ - α₂)))

        snum, rnum = splitsquare(s1n*s2n)
        sden, rden = splitsquare(s1d)
        snum, sden = divgcd!(snum, sden)
        rnum, rden = divgcd!(rnum, rden)
        s = Base.unsafe_rational(convert(BigInt, snum), convert(BigInt, sden))
        r = Base.unsafe_rational(convert(BigInt, rnum), convert(BigInt, rden))
        s *= compute3jseries(β₁, β₂, β₃, α₁, α₂)
        Wigner3j[(β₁, β₂, β₃, α₁, α₂)] = (r,s)
    end
    return _convert(T, sgn*s)*convert(T, signedroot(r))
end

"""
    clebschgordan(T::Type{<:Real} = RationalRoot{BigInt}, j₁, m₁, j₂, m₂, j₃, m₃ = m₁+m₂) -> ::T

Compute the value of the Clebsch-Gordan coefficient <j₁, m₁; j₂, m₂ | j₃, m₃ >
as a type `T` real number. By default, `T = RationalRoot{BigInt}` and `m₃ = m₁+m₂`.

Returns `zero(T)` if the triangle condition `δ(j₁, j₂, j₃)` is not satisfied, but
throws a `DomainError` if the `jᵢ`s and `mᵢ`s are not (half)integer or `abs(mᵢ) > jᵢ`.
"""
clebschgordan(j₁, m₁, j₂, m₂, j₃, m₃ = m₁+m₂) =
    clebschgordan(RRBig, j₁, m₁, j₂, m₂, j₃, m₃)
function clebschgordan(T::Type{<:Real}, j₁, m₁, j₂, m₂, j₃, m₃ = m₁+m₂)
    s = wigner3j(T, j₁, j₂, j₃, m₁, m₂, -m₃)
    iszero(s) && return s
    s *= convert(T, signedroot(RRBig, j₃+j₃+one(j₃)))
    return isodd(convert(Int,j₁ - j₂ + m₃)) ? -s : s
end

"""
    racahV(T::Type{<:Real} = RationalRoot{BigInt}, j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m₂) -> ::T

Compute the value of Racah's V-symbol
V(j₁, j₂, j₃; m₁, m₂, m₃) = (-1)^(-j₁+j₂+j₃) * ⎛ j₁  j₂  j₃ ⎞
                                               ⎝ m₁  m₂  m₃ ⎠
as a type `T` real number. By default, `T = RationalRoot{BigInt}` and `m₃ = -m₁-m₂`.

Returns `zero(T)` if the triangle condition `δ(j₁, j₂, j₃)` is not satisfied, but
throws a `DomainError` if the `jᵢ`s and `mᵢ`s are not (half)integer or `abs(mᵢ) > jᵢ`.
"""
racahV(j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m₂) = racahV(RRBig, j₁, j₂, j₃, m₁, m₂, m₃)
function racahV(T::Type{<:Real}, j₁, j₂, j₃, m₁, m₂, m₃ = -m₁-m₂)
    s = wigner3j(T, j₁, j₂, j₃, m₁, m₂, m₃)
    return isodd(convert(Int, -j₁ + j₂ + j₃)) ? -s : s
end

"""
    wigner6j(T::Type{<:Real} = RationalRoot{BigInt}, j₁, j₂, j₃, j₄, j₅, j₆) -> ::T

Compute the value of the Wigner-6j symbol
    _⎧ j₁  j₂  j₃ ⎫_
     ⎩ j₄  j₅  j₆ ⎭
as a type `T` real number. By default, `T = RationalRoot{BigInt}`.

Returns `zero(T)` if any of triangle conditions `δ(j₁, j₂, j₃)`, `δ(j₁, j₆, j₅)`,
`δ(j₂, j₄, j₆)`, `δ(j₃, j₄, j₅)` are not satisfied, but throws a `DomainError` if
the `jᵢ`s are not integer or halfinteger.
"""
wigner6j(j₁, j₂, j₃, j₄, j₅, j₆) = wigner6j(RRBig, j₁, j₂, j₃, j₄, j₅, j₆)
function wigner6j(T::Type{<:Real}, j₁, j₂, j₃, j₄, j₅, j₆)
    # check validity of `jᵢ`s
    for jᵢ in (j₁, j₂, j₃, j₄, j₅, j₆)
        (ishalfinteger(jᵢ) && jᵢ >= zero(jᵢ)) || throw(DomainError("invalid jᵢ", jᵢ))
    end
    return _wigner6j(T, HalfInteger.((j₁, j₂, j₃, j₄, j₅, j₆))...)
end

function _wigner6j(T::Type{<:Real}, j₁::HalfInteger, j₂::HalfInteger, j₃::HalfInteger,
                                    j₄::HalfInteger, j₅::HalfInteger, j₆::HalfInteger)
    α̂₁ = (j₁, j₂, j₃)
    α̂₂ = (j₁, j₆, j₅)
    α̂₃ = (j₂, j₄, j₆)
    α̂₄ = (j₃, j₄, j₅)

    # check triangle conditions
    if !(δ(α̂₁...) && δ(α̂₂...) && δ(α̂₃...) && δ(α̂₄...))
        return zero(T)
    end
    # reduce
    α₁ = convert(UInt, +(α̂₁...))
    α₂ = convert(UInt, +(α̂₂...))
    α₃ = convert(UInt, +(α̂₃...))
    α₄ = convert(UInt, +(α̂₄...))
    β₁ = convert(UInt, j₁+j₂+j₄+j₅)
    β₂ = convert(UInt, j₁+j₃+j₄+j₆)
    β₃ = convert(UInt, j₂+j₃+j₅+j₆)

    # we should have αᵢ < βⱼ, ∀ i, j and ∑ᵢ αᵢ = ∑ⱼ βⱼ
    # now order them as β₁ >= β₂ >= β₃ >= α₁ >= α₂ >= α₃ >= α₄
    (β₁, β₂, β₃, α₁, α₂, α₃, α₄) = reorder6j(β₁, β₂, β₃, α₁, α₂, α₃, α₄)

    # dictionary lookup or compute
    if haskey(Wigner6j, (β₁, β₂, β₃, α₁, α₂, α₃))
        r, s = Wigner6j[(β₁, β₂, β₃, α₁, α₂, α₃)]
    else
        # order irrelevant: product remains the same under action of reorder6j
        n₁, d₁ = Δ²(α̂₁...)
        n₂, d₂ = Δ²(α̂₂...)
        n₃, d₃ = Δ²(α̂₃...)
        n₄, d₄ = Δ²(α̂₄...)

        snum, rnum = splitsquare(n₁ * n₂ * n₃ * n₄)
        sden, rden = splitsquare(d₁ * d₂ * d₃ * d₄)
        snum, sden = divgcd!(snum, sden)
        rnum, rden = divgcd!(rnum, rden)
        s = Base.unsafe_rational(convert(BigInt, snum), convert(BigInt, sden))
        r = Base.unsafe_rational(convert(BigInt, rnum), convert(BigInt, rden))
        s *= compute6jseries(β₁, β₂, β₃, α₁, α₂, α₃, α₄)

        Wigner6j[(β₁, β₂, β₃, α₁, α₂, α₃)] = (r, s)
    end
    return _convert(T, s) * convert(T, signedroot(r))
end

"""
    racahW(T::Type{<:Real} = RationalRoot{BigInt}, j₁, j₂, J, j₃, J₁₂, J₂₃) -> ::T

Compute the value of Racah's W coefficient
`W(j₁, j₂, J, j₃; J₁₂, J₂₃) = <(j₁,(j₂j₃)J₂₃)J | ((j₁j₂)J₁₂,j₃)J> / sqrt((2J₁₂+1)*(2J₁₃+1))`
as a type `T` real number. By default, `T = RationalRoot{BigInt}`.

Returns `zero(T)` if any of triangle conditions `δ(j₁, j₂, J₁₂)`, `δ(j₂, j₃, J₂₃)`,
`δ(j₁, J₂₃, J)`, `δ(J₁₂, j₃, J)` are not satisfied, but throws a `DomainError` if
the `jᵢ`s and `J`s are not integer or halfinteger.
"""
racahW(j₁, j₂, J, j₃, J₁₂, J₂₃) = racahW(RRBig, j₁, j₂, J, j₃, J₁₂, J₂₃)
function racahW(T::Type{<:Real}, j₁, j₂, J, j₃, J₁₂, J₂₃)
    s = wigner6j(T, j₁, j₂, J₁₂, j₃, J, J₂₃)
    if !iszero(s) && isodd(convert(Int, j₁ + j₂ + j₃ + J))
        return -s
    else
        return s
    end
end

# COMPUTATIONAL ROUTINES
#------------------------
# squared triangle coefficient
function Δ²(j₁, j₂, j₃)
    # also checks the triangle conditions by converting to unsigned integer:
    n1 = copy(primefactorial( convert(UInt, + j₁ + j₂ - j₃) ))
    n2 = primefactorial( convert(UInt, + j₁ - j₂ + j₃) )
    n3 = primefactorial( convert(UInt, - j₁ + j₂ + j₃) )
    num = mul!(mul!(n1, n2), n3)
    den = copy(primefactorial( convert(UInt, j₁ + j₂ + j₃ + 1) ))
    # result
    return divgcd!(num, den)
end

# reorder parameters determining the 3j symbol to canonical order:
# j₁ >= j₂ >= j₃ and m₁ >= 0 or m₁ == 0 && m₂ >= 0
function reorder3j(j₁, j₂, j₃, m₁, m₂, m₃, sign = one(Int8))
    if j₁ < j₂
        return reorder3j(j₂, j₁, j₃, m₂, m₁, m₃, -sign)
    elseif j₂ < j₃
        return reorder3j(j₁, j₃, j₂, m₁, m₃, m₂, -sign)
    elseif m₁ < zero(m₁)
        return reorder3j(j₁, j₂, j₃, -m₁, -m₂, -m₃, -sign)
    elseif iszero(m₁) && m₂ < zero(m₂)
        return reorder3j(j₁, j₂, j₃, -m₁, -m₂, -m₃, -sign)
    else
        # sign doesn't matter if total J=j₁ + j₂ + j₃ is even
        if iseven(convert(UInt,j₁ + j₂ + j₃))
            sign = one(sign)
        end
        return (j₁, j₂, j₃, m₁, m₂, m₃, sign)
    end
end

# reorder parameters determining the 6j symbol to canonical order:
# β₁ >= β₂ >= β₃ >= α₁ >= α₂ >= α₃ >= α₄
function reorder6j(β₁, β₂, β₃, α₁, α₂, α₃, α₄)
    if β₁ < β₂
        return reorder6j(β₂, β₁, β₃, α₁, α₂, α₃, α₄)
    elseif β₂ < β₃
        return reorder6j(β₁, β₃, β₂, α₁, α₂, α₃, α₄)
    elseif α₁ < α₂
        return reorder6j(β₁, β₂, β₃, α₂, α₁, α₃, α₄)
    elseif α₂ < α₃
        return reorder6j(β₁, β₂, β₃, α₁, α₃, α₂, α₄)
    elseif α₃ < α₄
        return reorder6j(β₁, β₂, β₃, α₁, α₂, α₄, α₃)
    else
        return (β₁, β₂, β₃, α₁, α₂, α₃, α₄)
    end
end

# compute the sum appearing in the 3j symbol
function compute3jseries(β₁, β₂, β₃, α₁, α₂)
    krange = max(α₁, α₂, zero(α₁)):min(β₁, β₂, β₃)
    T = PrimeFactorization{eltype(eltype(factorialtable))}

    nums = Vector{T}(undef, length(krange))
    dens = Vector{T}(undef, length(krange))
    for (i, k) in enumerate(krange)
        num = iseven(k) ? one(T) : -one(T)
        den = copy(primefactorial(k))
        den = mul!(mul!(den, primefactorial(k-α₁)), primefactorial(k-α₂))
        den = mul!(mul!(mul!(den, primefactorial(β₁-k)),
                                    primefactorial(β₂-k)),
                                        primefactorial(β₃-k))
        nums[i], dens[i] = num, den
    end
    den = commondenominator!(nums, dens)
    totalnum = sumlist!(nums)
    totalden = convert(BigInt, den)
    for n = 1:length(den.powers)
        p = bigprime(n)
        while den.powers[n] > 0
            q, r = divrem(totalnum, p)
            if iszero(r)
                totalnum = q
                den.powers[n] -= 1
            else
                break
            end
        end
    end
    totalden = convert(BigInt, den)
    return Base.unsafe_rational(totalnum, totalden)
end

# compute the sum appearing in the 6j symbol
function compute6jseries(β₁, β₂, β₃, α₁, α₂, α₃, α₄)
    krange = max(α₁, α₂, α₃, α₄):min(β₁, β₂, β₃)
    T = PrimeFactorization{eltype(eltype(factorialtable))}

    nums = Vector{T}(undef, length(krange))
    dens = Vector{T}(undef, length(krange))
    for (i, k) in enumerate(krange)
        num = iseven(k) ? copy(primefactorial(k+1)) : neg!(copy(primefactorial(k+1)))
        den = copy(primefactorial(k-α₁))
        den = mul!(mul!(mul!(den, primefactorial(k-α₂)),
                                    primefactorial(k-α₃)),
                                        primefactorial(k-α₄))
        den = mul!(mul!(mul!(den, primefactorial(β₁-k)),
                                    primefactorial(β₂-k)),
                                        primefactorial(β₃-k))
        nums[i], dens[i] = divgcd!(num, den)
    end
    den = commondenominator!(nums, dens)
    totalnum = sumlist!(nums)
    for n = 1:length(den.powers)
        p = bigprime(n)
        while den.powers[n] > 0
            q, r = divrem(totalnum, p)
            if iszero(r)
                totalnum = q
                den.powers[n] -= 1
            else
                break
            end
        end
    end
    totalden = convert(BigInt, den)
    return Base.unsafe_rational(totalnum, totalden)
end

function _precompile_()
    @assert precompile(wigner3j, (Type{Float64}, Int, Int, Int, Int, Int, Int))
    @assert precompile(wigner6j, (Type{Float64}, Int, Int, Int, Int, Int, Int))
    @assert precompile(wigner3j, (Type{BigFloat}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}))
    @assert precompile(wigner6j, (Type{BigFloat}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}))
    @assert precompile(wigner3j, (Type{RRBig}, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt))
    @assert precompile(wigner6j, (Type{RRBig}, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt))
end
_precompile_()

end # module


text:
#! /usr/bin/julia

# Rosetta Code, Cholesky decomposition

a = [25 15 5; 15 18 0; -5 0 11]
b = [18 22 54 22; 22 70 86 62; 54 86 174 134; 42 62 134 106]

println(a, "\n => \n", chol(a, :L))
println(b, "\n => \n", chol(b, :L))

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using JSON, Mustache

"""
    convert_doc(infile::AbstractString, outfile::AbstractString; outformat::Union{Nothing,AbstractString} = nothing)

Convert Weave documents between different formats

- `infile`: Path of the input document
- `outfile`: Path of the output document
- `outformat = nothing`: Output document format (optional). By default (i.e. given `nothing`) Weave will try to automatically detect it from the `outfile`'s extension. You can also specify either of `"script"`, `"markdown"`, `"notebook"`, or `"noweb"`
"""
function convert_doc(
    infile::AbstractString,
    outfile::AbstractString;
    outformat::Union{Nothing,AbstractString} = nothing,
)
    doc = WeaveDoc(infile)

    if isnothing(outformat)
        ext = lowercase(splitext(outfile)[2])
        outformat =
            ext == ".jl" ? "script" :
            ext == ".jmd" ?  "markdown" :
            ext == ".ipynb" ? "notebook" :
            "noweb" # fallback
    end

    converted = _convert_doc(doc, outformat)

    open(outfile, "w") do f
        write(f, converted)
    end
    return outfile
end

function _convert_doc(doc, outformat)
    outformat == "script" ? convert_to_script(doc) :
    outformat == "markdown" ? convert_to_markdown(doc) :
    outformat == "notebook" ? convert_to_notebook(doc) :
    convert_to_noweb(doc)
end

function convert_to_script(doc)
    output = ""
    for chunk in doc.chunks
        if typeof(chunk) == Weave.DocChunk
            content = join([repr(c) for c in chunk.content], "")
            output *= join(["#' " * s for s in split(content, "\n")], "\n")
        else
            output *= "\n#+ "
            isempty(chunk.optionstring) || (output *= strip(chunk.optionstring))
            output *= "\n\n" * lstrip(chunk.content)
            output *= "\n"
        end
    end

    return output
end

function convert_to_markdown(doc)
    output = ""
    for chunk in doc.chunks
        if isa(chunk, DocChunk)
            output *= join([repr(c) for c in chunk.content], "")
        else
            output *= "\n" * "```julia"
            isempty(chunk.optionstring) || (output *= ";" * chunk.optionstring)
            output *= "\n" * lstrip(chunk.content)
            output *= "```\n"
        end
    end

    return output
end

function convert_to_notebook(doc)
    nb = Dict()
    nb["nbformat"] = 4
    nb["nbformat_minor"] = 2
    metadata = Dict()
    kernelspec = Dict()
    kernelspec["language"] = "julia"
    kernelspec["name"] = "julia-$(VERSION.major).$(VERSION.minor)"
    kernelspec["display_name"] = "Julia $(VERSION.major).$(VERSION.minor).$(VERSION.patch)"
    metadata["kernelspec"] = kernelspec
    language_info = Dict()
    language_info["file_extension"] = ".jl"
    language_info["mimetype"] = "application/julia"
    language_info["name"] = "julia"
    language_info["version"] = "$(VERSION.major).$(VERSION.minor).$(VERSION.patch)"
    metadata["language_info"] = language_info
    cells = []
    ex_count = 1

    for chunk in doc.chunks
        if isa(chunk, DocChunk)
            push!(
                cells,
                Dict(
                    "cell_type" => "markdown",
                    "metadata" => Dict(),
                    "source" => [strip(join([repr(c) for c in chunk.content], ""))],
                ),
            )
        else
            push!(
                cells,
                Dict(
                    "cell_type" => "code",
                    "metadata" => Dict(),
                    "source" => [strip(chunk.content)],
                    "execution_count" => nothing,
                    "outputs" => [],
                ),
            )
        end
    end

    nb["cells"] = cells
    nb["metadata"] = metadata

    json_nb = JSON.json(nb, 2)
    return json_nb
end

function convert_to_noweb(doc)
    output = ""
    for chunk in doc.chunks
        if isa(chunk, DocChunk)
            output *= join([repr(c) for c in chunk.content], "")
        else
            output *= "\n" * "<<"
            isempty(chunk.optionstring) || (output *= strip(chunk.optionstring))
            output *= ">>="
            output *= "\n" * lstrip(chunk.content)
            output *= "@\n"
        end
    end

    return output
end

Base.repr(c::InlineText) = c.content
Base.repr(c::InlineCode) = "`j $(c.content)`"


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## shell-like command parsing ##

const shell_special = "#{}()[]<>|&*?~;"

# strips the end but respects the space when the string ends with "\\ "
function rstrip_shell(s::AbstractString)
    c_old = nothing
    for (i, c) in Iterators.reverse(pairs(s))
        i::Int; c::AbstractChar
        ((c == '\\') && c_old == ' ') && return SubString(s, 1, i+1)
        isspace(c) || return SubString(s, 1, i)
        c_old = c
    end
    SubString(s, 1, 0)
end

function shell_parse(str::AbstractString, interpolate::Bool=true;
                     special::AbstractString="", filename="none")
    s = SubString(str, firstindex(str))
    s = rstrip_shell(lstrip(s))

    # N.B.: This is used by REPLCompletions
    last_parse = 0:-1
    isempty(s) && return interpolate ? (Expr(:tuple,:()),last_parse) : ([],last_parse)

    in_single_quotes = false
    in_double_quotes = false

    args = []
    arg = []
    i = firstindex(s)
    st = Iterators.Stateful(pairs(s))

    function push_nonempty!(list, x)
        if !isa(x,AbstractString) || !isempty(x)
            push!(list, x)
        end
        return nothing
    end
    function consume_upto!(list, s, i, j)
        push_nonempty!(list, s[i:prevind(s, j)::Int])
        something(peek(st), lastindex(s)::Int+1 => '\0').first::Int
    end
    function append_2to1!(list, innerlist)
        if isempty(innerlist); push!(innerlist, ""); end
        push!(list, copy(innerlist))
        empty!(innerlist)
    end

    C = eltype(str)
    P = Pair{Int,C}
    for (j, c) in st
        j, c = j::Int, c::C
        if !in_single_quotes && !in_double_quotes && isspace(c)
            i = consume_upto!(arg, s, i, j)
            append_2to1!(args, arg)
            while !isempty(st)
                # We've made sure above that we don't end in whitespace,
                # so updating `i` here is ok
                (i, c) = peek(st)::P
                isspace(c) || break
                popfirst!(st)
            end
        elseif interpolate && !in_single_quotes && c == '$'
            i = consume_upto!(arg, s, i, j)
            isempty(st) && error("\$ right before end of command")
            stpos, c = popfirst!(st)::P
            isspace(c) && error("space not allowed right after \$")
            if startswith(SubString(s, stpos), "var\"")
                # Disallow var"#" syntax in cmd interpolations.
                # TODO: Allow only identifiers after the $ for consistency with
                # string interpolation syntax (see #3150)
                ex, j = :var, stpos+3
            else
                # use parseatom instead of parse to respect filename (#28188)
                ex, j = Meta.parseatom(s, stpos, filename=filename)
            end
            last_parse = (stpos:prevind(s, j)) .+ s.offset
            push_nonempty!(arg, ex)
            s = SubString(s, j)
            Iterators.reset!(st, pairs(s))
            i = firstindex(s)
        else
            if !in_double_quotes && c == '\''
                in_single_quotes = !in_single_quotes
                i = consume_upto!(arg, s, i, j)
            elseif !in_single_quotes && c == '"'
                in_double_quotes = !in_double_quotes
                i = consume_upto!(arg, s, i, j)
            elseif !in_single_quotes && c == '\\'
                if !isempty(st) && (peek(st)::P)[2] in ('\n', '\r')
                    i = consume_upto!(arg, s, i, j) + 1
                    if popfirst!(st)[2] == '\r' && (peek(st)::P)[2] == '\n'
                        i += 1
                        popfirst!(st)
                    end
                    while !isempty(st) && (peek(st)::P)[2] in (' ', '\t')
                        i = nextind(str, i)
                        _ = popfirst!(st)
                    end
                elseif in_double_quotes
                    isempty(st) && error("unterminated double quote")
                    k, c′ = peek(st)::P
                    if c′ == '"' || c′ == '$' || c′ == '\\'
                        i = consume_upto!(arg, s, i, j)
                        _ = popfirst!(st)
                    end
                else
                    isempty(st) && error("dangling backslash")
                    i = consume_upto!(arg, s, i, j)
                    _ = popfirst!(st)
                end
            elseif !in_single_quotes && !in_double_quotes && c in special
                error("parsing command `$str`: special characters \"$special\" must be quoted in commands")
            end
        end
    end

    if in_single_quotes; error("unterminated single quote"); end
    if in_double_quotes; error("unterminated double quote"); end

    push_nonempty!(arg, s[i:end])
    append_2to1!(args, arg)

    interpolate || return args, last_parse

    # construct an expression
    ex = Expr(:tuple)
    for arg in args
        push!(ex.args, Expr(:tuple, arg...))
    end
    return ex, last_parse
end

function shell_split(s::AbstractString)
    parsed = shell_parse(s, false)[1]
    args = String[]
    for arg in parsed
        push!(args, string(arg...))
    end
    args
end

function print_shell_word(io::IO, word::AbstractString, special::AbstractString = "")
    has_single = false
    has_special = false
    for c in word
        if isspace(c) || c=='\\' || c=='\'' || c=='"' || c=='$' || c in special
            has_special = true
            if c == '\''
                has_single = true
            end
        end
    end
    if isempty(word)
        print(io, "''")
    elseif !has_special
        print(io, word)
    elseif !has_single
        print(io, '\'', word, '\'')
    else
        print(io, '"')
        for c in word
            if c == '"' || c == '$'
                print(io, '\\')
            end
            print(io, c)
        end
        print(io, '"')
    end
    nothing
end

function print_shell_escaped(io::IO, cmd::AbstractString, args::AbstractString...;
                             special::AbstractString="")
    print_shell_word(io, cmd, special)
    for arg in args
        print(io, ' ')
        print_shell_word(io, arg, special)
    end
end
print_shell_escaped(io::IO; special::String="") = nothing

"""
    shell_escape(args::Union{Cmd,AbstractString...}; special::AbstractString="")

The unexported `shell_escape` function is the inverse of the unexported `shell_split` function:
it takes a string or command object and escapes any special characters in such a way that calling
`shell_split` on it would give back the array of words in the original command. The `special`
keyword argument controls what characters in addition to whitespace, backslashes, quotes and
dollar signs are considered to be special (default: none).

# Examples
```jldoctest
julia> Base.shell_escape("cat", "/foo/bar baz", "&&", "echo", "done")
"cat '/foo/bar baz' && echo done"

julia> Base.shell_escape("echo", "this", "&&", "that")
"echo this && that"
```
"""
shell_escape(args::AbstractString...; special::AbstractString="") =
    sprint((io, args...) -> print_shell_escaped(io, args..., special=special), args...)

function print_shell_escaped_posixly(io::IO, args::AbstractString...)
    first = true
    for arg in args
        first || print(io, ' ')
        # avoid printing quotes around simple enough strings
        # that any (reasonable) shell will definitely never consider them to be special
        have_single = false
        have_double = false
        function isword(c::AbstractChar)
            if '0' <= c <= '9' || 'a' <= c <= 'z' || 'A' <= c <= 'Z'
                # word characters
            elseif c == '_' || c == '/' || c == '+' || c == '-'
                # other common characters
            elseif c == '\''
                have_single = true
            elseif c == '"'
                have_double && return false # switch to single quoting
                have_double = true
            elseif !first && c == '='
                # equals is special if it is first (e.g. `env=val ./cmd`)
            else
                # anything else
                return false
            end
            return true
        end
        if isempty(arg)
            print(io, "''")
        elseif all(isword, arg)
            have_single && (arg = replace(arg, '\'' => "\\'"))
            have_double && (arg = replace(arg, '"' => "\\\""))
            print(io, arg)
        else
            print(io, '\'', replace(arg, '\'' => "'\\''"), '\'')
        end
        first = false
    end
end

"""
    shell_escape_posixly(args::Union{Cmd,AbstractString...})

The unexported `shell_escape_posixly` function
takes a string or command object and escapes any special characters in such a way that
it is safe to pass it as an argument to a posix shell.

# Examples
```jldoctest
julia> Base.shell_escape_posixly("cat", "/foo/bar baz", "&&", "echo", "done")
"cat '/foo/bar baz' '&&' echo done"

julia> Base.shell_escape_posixly("echo", "this", "&&", "that")
"echo this '&&' that"
```
"""
shell_escape_posixly(args::AbstractString...) =
    sprint(print_shell_escaped_posixly, args...)

"""
    shell_escape_csh(args::Union{Cmd,AbstractString...})
    shell_escape_csh(io::IO, args::Union{Cmd,AbstractString...})

This function quotes any metacharacters in the string arguments such
that the string returned can be inserted into a command-line for
interpretation by the Unix C shell (csh, tcsh), where each string
argument will form one word.

In contrast to a POSIX shell, csh does not support the use of the
backslash as a general escape character in double-quoted strings.
Therefore, this function wraps strings that might contain
metacharacters in single quotes, except for parts that contain single
quotes, which it wraps in double quotes instead. It switches between
these types of quotes as needed. Linefeed characters are escaped with
a backslash.

This function should also work for a POSIX shell, except if the input
string contains a linefeed (`"\\n"`) character.

See also: [`shell_escape_posixly`](@ref)
"""
function shell_escape_csh(io::IO, args::AbstractString...)
    first = true
    for arg in args
        first || write(io, ' ')
        first = false
        i = 1
        while true
            for (r,e) = (r"^[A-Za-z0-9/\._-]+\z" => "",
                         r"^[^']*\z" => "'", r"^[^\$\`\"]*\z" => "\"",
                         r"^[^']+"  => "'", r"^[^\$\`\"]+"  => "\"")
                if ((m = match(r, SubString(arg, i))) !== nothing)
                    write(io, e)
                    write(io, replace(m.match, '\n' => "\\\n"))
                    write(io, e)
                    i += ncodeunits(m.match)
                    break
                end
            end
            i <= lastindex(arg) || break
        end
    end
end
shell_escape_csh(args::AbstractString...) =
    sprint(shell_escape_csh, args...;
           sizehint = sum(sizeof.(args)) + length(args) * 3)

"""
    shell_escape_wincmd(s::AbstractString)
    shell_escape_wincmd(io::IO, s::AbstractString)

The unexported `shell_escape_wincmd` function escapes Windows `cmd.exe` shell
meta characters. It escapes `()!^<>&|` by placing a `^` in front. An `@` is
only escaped at the start of the string. Pairs of `"` characters and the
strings they enclose are passed through unescaped. Any remaining `"` is escaped
with `^` to ensure that the number of unescaped `"` characters in the result
remains even.

Since `cmd.exe` substitutes variable references (like `%USER%`) _before_
processing the escape characters `^` and `"`, this function makes no attempt to
escape the percent sign (`%`), the presence of `%` in the input may cause
severe breakage, depending on where the result is used.

Input strings with ASCII control characters that cannot be escaped (NUL, CR,
LF) will cause an `ArgumentError` exception.

The result is safe to pass as an argument to a command call being processed by
`CMD.exe /S /C " ... "` (with surrounding double-quote pair) and will be
received verbatim by the target application if the input does not contain `%`
(else this function will fail with an ArgumentError). The presence of `%` in
the input string may result in command injection vulnerabilities and may
invalidate any claim of suitability of the output of this function for use as
an argument to cmd (due to the ordering described above), so use caution when
assembling a string from various sources.

This function may be useful in concert with the `windows_verbatim` flag to
[`Cmd`](@ref) when constructing process pipelines.

```julia
wincmd(c::String) =
   run(Cmd(Cmd(["cmd.exe", "/s /c \\" \$c \\""]);
           windows_verbatim=true))
wincmd_echo(s::String) =
   wincmd("echo " * Base.shell_escape_wincmd(s))
wincmd_echo("hello \$(ENV["USERNAME"]) & the \\"whole\\" world! (=^I^=)")
```

But take note that if the input string `s` contains a `%`, the argument list
and echo'ed text may get corrupted, resulting in arbitrary command execution.
The argument can alternatively be passed as an environment variable, which
avoids the problem with `%` and the need for the `windows_verbatim` flag:

```julia
cmdargs = Base.shell_escape_wincmd("Passing args with %cmdargs% works 100%!")
run(setenv(`cmd /C echo %cmdargs%`, "cmdargs" => cmdargs))
```

!warning
    The argument parsing done by CMD when calling batch files (either inside
    `.bat` files or as arguments to them) is not fully compatible with the
    output of this function. In particular, the processing of `%` is different.

!important
    Due to a peculiar behavior of the CMD parser/interpreter, each command
    after a literal `|` character (indicating a command pipeline) must have
    `shell_escape_wincmd` applied twice since it will be parsed twice by CMD.
    This implies ENV variables would also be expanded twice!
    For example:
    ```julia
    to_print = "All for 1 & 1 for all!"
    to_print_esc = Base.shell_escape_wincmd(Base.shell_escape_wincmd(to_print))
    run(Cmd(Cmd(["cmd", "/S /C \\" break | echo \$(to_print_esc) \\""]), windows_verbatim=true))
    ```

With an I/O stream parameter `io`, the result will be written there,
rather than returned as a string.

See also [`escape_microsoft_c_args`](@ref), [`shell_escape_posixly`](@ref).

# Example
```jldoctest
julia> Base.shell_escape_wincmd("a^\\"^o\\"^u\\"")
"a^^\\"^o\\"^^u^\\""
```
"""
function shell_escape_wincmd(io::IO, s::AbstractString)
    # https://stackoverflow.com/a/4095133/1990689
    occursin(r"[\r\n\0]", s) &&
        throw(ArgumentError("control character unsupported by CMD.EXE"))
    i = 1
    len = ncodeunits(s)
    if len > 0 && s[1] == '@'
        write(io, '^')
    end
    while i <= len
        c = s[i]
        if c == '"' && (j = findnext('"', s, nextind(s,i))) !== nothing
            write(io, SubString(s,i,j))
            i = j
        else
            if c in ('"', '(', ')', '!', '^', '<', '>', '&', '|')
                write(io, '^', c)
            else
                write(io, c)
            end
        end
        i = nextind(s,i)
    end
end
shell_escape_wincmd(s::AbstractString) = sprint(shell_escape_wincmd, s;
                                                sizehint = 2*sizeof(s))

"""
    escape_microsoft_c_args(args::Union{Cmd,AbstractString...})
    escape_microsoft_c_args(io::IO, args::Union{Cmd,AbstractString...})

Convert a collection of string arguments into a string that can be
passed to many Windows command-line applications.

Microsoft Windows passes the entire command line as a single string to
the application (unlike POSIX systems, where the shell splits the
command line into a list of arguments). Many Windows API applications
(including julia.exe), use the conventions of the [Microsoft C/C++
runtime](https://docs.microsoft.com/en-us/cpp/c-language/parsing-c-command-line-arguments)
to split that command line into a list of strings.

This function implements an inverse for a parser compatible with these rules.
It joins command-line arguments to be passed to a Windows
C/C++/Julia application into a command line, escaping or quoting the
meta characters space, TAB, double quote and backslash where needed.

See also [`shell_escape_wincmd`](@ref), [`escape_raw_string`](@ref).
"""
function escape_microsoft_c_args(io::IO, args::AbstractString...)
    # http://daviddeley.com/autohotkey/parameters/parameters.htm#WINCRULES
    first = true
    for arg in args
        if first
            first = false
        else
            write(io, ' ')  # separator
        end
        if isempty(arg) || occursin(r"[ \t\"]", arg)
            # Julia raw strings happen to use the same escaping convention
            # as the argv[] parser in Microsoft's C runtime library.
            write(io, '"')
            escape_raw_string(io, arg)
            write(io, '"')
        else
            write(io, arg)
        end
    end
end
escape_microsoft_c_args(args::AbstractString...) =
    sprint(escape_microsoft_c_args, args...;
           sizehint = (sum(sizeof.(args)) + 3*length(args)))

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
module Rsvg

include("../deps/deps.jl")

using Cairo
using Gtk

include("types.jl")
include("calls.jl")

export RsvgDimensionData, RsvgHandle

Rsvg.set_default_dpi(72.0)

end # module
#


text:
# Modified from https://github.com/r3tex/ObjectDetector.jl/blob/master/dev/artifacts/generate_artifacts.jl
# First, run `download_data.jl` from the `/dev` directory to download the files into `/data`
# Then run this file to upload them to Github releases and populate `/Artifacts.toml`

using Pkg.Artifacts, LibGit2, ghr_jll

deploy = true
if deploy && !haskey(ENV, "GITHUB_TOKEN")
    error("For automatic github deployment, export GITHUB_TOKEN!")
end

data_dir = joinpath(@__DIR__, "..", "data")

if deploy
    # Where we will put our tarballs
    tmpdir = mktempdir()

    function get_git_remote_url(repo_path::String)
        repo = LibGit2.GitRepo(repo_path)
        origin = LibGit2.get(LibGit2.GitRemote, repo, "origin")
        return LibGit2.url(origin)
    end

    # Try to detect where we should upload these weights to (or just override
    # as shown in the commented-out line)
    origin_url = get_git_remote_url(dirname(@__DIR__))
    deploy_repo = "$(basename(dirname(origin_url)))/$(basename(origin_url))"
    #deploy_repo = "staticfloat/ObjectDetector.jl"

    tag = "traineddata"
end

for dir in readdir(data_dir)
    weightsfiles = filter(x -> endswith(x, ".traineddata"),
                          readdir(joinpath(data_dir, dir)))

    @info("Generating artifact for $(dir)")
    # Create a local artifact
    hash = create_artifact() do artifact_dir
        # Copy in weights
        for weightsfile in weightsfiles
            cp(joinpath(data_dir, dir, weightsfile), joinpath(artifact_dir, weightsfile))
        end
        # Copy in the license
        cp(joinpath(data_dir, dir, "LICENSE"), joinpath(artifact_dir, "LICENSE"))
        # Copy in the `pdf.ttf` file
        return cp(joinpath(data_dir, dir, "pdf.ttf"), joinpath(artifact_dir, "pdf.ttf"))
    end

    # Spit tarballs to be hosted out to local temporary directory:
    if deploy
        tarball_hash = archive_artifact(hash, joinpath(tmpdir, "$(dir).tar.gz"))

        # Calculate tarball url
        tarball_url = "https://github.com/$(deploy_repo)/releases/download/$(tag)/$(dir).tar.gz"

        # Bind this to an Artifacts.toml file
        @info("Binding $(dir) in Artifacts.toml...")
        bind_artifact!(joinpath(@__DIR__, "..", "Artifacts.toml"), dir, hash;
                       download_info=[(tarball_url, tarball_hash)], lazy=true, force=true)
    end
end

if deploy
    # Upload tarballs to a special github release
    @info("Uploading tarballs to $(deploy_repo) tag `$(tag)`")
    ghr() do ghr_exe
        return run(`$ghr_exe -replace -u $(dirname(deploy_repo)) -r $(basename(deploy_repo)) $(tag) $(tmpdir)`)
    end

    @info("Artifacts.toml file now contains all bound artifact names")
    rm(tmpdir; force=true, recursive=true)
end


text:
"like `copy!` but allowing for heterogeneous matrix types"
inject!(d,s) = copy!(d,s)               # fallback method

function inject!(d::UpperTriangular,s::UpperTriangular)
    (n = size(s,2)) == size(d,2) || throw(DimensionMismatch(""))
    for j in 1:n
        inject!(sub(d,1:j,j),sub(s,1:j,j))
    end
    d
end

function inject!{T<:Real}(d::AbstractMatrix{T}, s::Diagonal{T})
    sd = s.diag
    length(sd) == size(d,2) || throw(DimensionMismatch(""))
    fill!(d,0.)
    @inbounds for i in eachindex(sd)
        d[i,i] = sd[i]
    end
    d
end

inject!(d::Diagonal{Float64},s::Diagonal{Float64}) = (copy!(d.diag,s.diag);d)

function inject!(d::SparseMatrixCSC{Float64},s::SparseMatrixCSC{Float64})
    m,n = size(d)
    size(d) == size(s) || throw(DimensionMismatch(""))
    if nnz(d) == nnz(s)
        copy!(nonzeros(d),nonzeros(s))
        return d
    end
    drv = rowvals(d); srv = rowvals(s); dnz = nonzeros(d); snz = nonzeros(s)
    fill!(dnz,0.)
    for j in 1:n
        dnzr = nzrange(d,j)
        dnzrv = sub(drv,dnzr)
        snzr = nzrange(s,j)
        if length(snzr) == length(dnzr) && all(dnzrv .== sub(srv,snzr))
            copy!(sub(dnz,dnzr),sub(snz,snzr))
        else
            for k in snzr
                ssr = srv[k]
                kk = searchsortedfirst(dnzrv,ssr)
                kk > length(dnzrv) || dnzrv[kk] != ssr || error("cannot inject sparse s into sparse d")
                dnz[dnzr[kk]] = snz[k]
            end
        end
    end
    d
end


text:
"""
    my_f(x,y)

Returns double the number 'x' plus triple the number 'y'.

"""
function my_f(x, y)
    return 2x + 3y
end

"""
    my_f_add(x,y)

Returns double the number 'x' plus double the number 'y'.

"""
function my_f_add(x, y)
    return 2x + 2y
end

"""
    add_exported(x, y)

Dummy function, exported

# Examples
```jldoctest
julia> add_exported(1, 2)
3
```
"""
function add_exported(x::Number, y::Number)
    return x + y
end

#"""
#    add_not_exported(x, y)

#Dummy function, not exported

## Examples
#```jldoctest
#julia> packageTrial.add_not_exported(1, 2)
#3
#```
#"""
#function add_not_exported(x::Number, y::Number)
#    return x + y
#end


text:
function construct_model(t_train, t_test, u0_2, u0_3; hidden_dim = 10, test_model = true)
    ann = Chain(Dense(4, hidden_dim, relu), Dense(hidden_dim, 1))
    p1, re = Flux.destructure(ann)
    p2 = Float64[0.5, 0.03]
    p3 = [p1; p2]

    ps = Flux.params(p3)

    u0 = Float64[25.0e6, u0_2, u0_3, 10]

    function QSIR(du, u, p, t)
        β = abs(p[length(p1) + 1])
        γ = 0.007627
        du[1]=  -β * u[1] * u[2] / u0[1]
        du[2] = β * u[1] * u[2] / u0[1] - γ * u[2] - abs(re(p[1:length(p1)])(u)[1]) * u[2] / u0[1]
        du[3] = γ * u[2]
        du[4] = abs(re(p[1:length(p1)])(u)[1] * u[2] / u0[1])
    end

    tspan = (1.0, maximum(t_test))
    prob = ODEProblem(QSIR, u0, tspan, p3)

    if test_model
        sol = Array(concrete_solve(prob, Rosenbrock23(autodiff = false), u0, p3, saveat=t_train))
    end

    return ps, re, u0, prob, p3, length(p1)
end