text:
function bucr2_model(day::Int64, hour::Int64, df_gens::DataFrame,
  df_peakers::DataFrame, df_storage::DataFrame, fuelprice::Matrix{Float64},
  fuelprice_peakers::Matrix{Float64}, map_gens::Matrix{Int64},
  map_peakers::Matrix{Int64}, map_storage::Matrix{Int64}, fucr::UC_Results,
  bucr1::UC_Results, bucr2::UC_Results, sucr_to_bucr2::UC_Results,
  bucr2_prepoc::DemandHrPreprocGens, bucr_commit_LB::Vector{Int64}, bucr_commit_UB::Vector{Int64},
  bucr_commit_peaker_LB::Vector{Int64}, bucr_commit_peaker_UB::Vector{Int64})

	t1_BUCR2model = time_ns()

	BUCR2model = direct_model(CPLEX.Optimizer())
	set_optimizer_attribute(BUCR2model, "CPX_PARAM_EPGAP", SOLVER_EPGAP)

	# Declaring the decision variables for conventional generators
	@variable(BUCR2model, bucrm_genOnOff[1:GENS], Bin) #Bin
	@variable(BUCR2model, bucrm_genStartUp[1:GENS], Bin) # startup variable
	@variable(BUCR2model, bucrm_genShutDown[1:GENS], Bin) # shutdown variable
	@variable(BUCR2model, bucrm_genOut[1:GENS]>=0) # Generator's output schedule
	# Generator's output schedule from each block (block refers to IHR curve blocks/segments)
	@variable(BUCR2model, bucrm_genOut_Block[1:GENS, 1:BLOCKS]>=0)
	@variable(BUCR2model, bucrm_TotGenVioP[g=1:GENS] >= 0)
	@variable(BUCR2model, bucrm_TotGenVioN[g=1:GENS] >= 0)
	@variable(BUCR2model, bucrm_MaxGenVioP[g=1:GENS] >= 0)
	@variable(BUCR2model, bucrm_MinGenVioP[g=1:GENS] >= 0)

	# Declaring the decision variables for peakers
	@variable(BUCR2model, 0<=bucrm_peakerOnOff[1:PEAKERS]<=1) #Bin, relaxed
	@variable(BUCR2model, 0<=bucrm_peakerStartUp[1:PEAKERS]<=1) # startup variable
	@variable(BUCR2model, 0<=bucrm_peakerShutDown[1:PEAKERS]<=1) # shutdown variable
	@variable(BUCR2model, bucrm_peakerOut[1:PEAKERS]>=0) # Generator's output schedule
	@variable(BUCR2model, bucrm_peakerOut_Block[1:PEAKERS, 1:BLOCKS]>=0) # Generator's output schedule from each block (block rfers to IHR curve blocks/segments)

	#@variable(BUCR2model, bucrm_genResUp[1:GENS]>=0) # Generators' up reserve schedule
	#@variable(BUCR2model, bucrm_genResDn[1:GENS]>=0) # Generator's down rserve schedule

	# declaring decision variables for storage Units
	@variable(BUCR2model, bucrm_storgChrg[1:STORG_UNITS], Bin) #Bin variable equal to 1 if unit runs in the charging mode
	@variable(BUCR2model, bucrm_storgDisc[1:STORG_UNITS], Bin) #Bin variable equal to 1 if unit runs in the discharging mode
	@variable(BUCR2model, bucrm_storgIdle[1:STORG_UNITS], Bin) ##Bin variable equal to 1 if unit runs in the idle mode
	@variable(BUCR2model, bucrm_storgChrgPwr[1:STORG_UNITS]>=0) #Chargung power
	@variable(BUCR2model, bucrm_storgDiscPwr[1:STORG_UNITS]>=0) # Discharging Power
	@variable(BUCR2model, bucrm_storgSOC[1:STORG_UNITS]>=0) # state of charge (stored energy level for storage unit at time t)
	#@variable(BUCR2model, bucrm_storgResUp[1:STORG_UNITS]>=0) # Scheduled up reserve
	#@variable(BUCR2model, bucrm_storgResDn[1:STORG_UNITS]>=0) # Scheduled down reserve

	# declaring decision variables for renewable generation
	@variable(BUCR2model, bucrm_solarG[1:N_ZONES]>=0) # solar energy schedules
	@variable(BUCR2model, bucrm_windG[1:N_ZONES]>=0) # wind energy schedules
	@variable(BUCR2model, bucrm_hydroG[1:N_ZONES]>=0) # hydro energy schedules
	@variable(BUCR2model, bucrm_solarGSpil[1:N_ZONES]>=0) # solar energy schedules
	@variable(BUCR2model, bucrm_windGSpil[1:N_ZONES]>=0) # wind energy schedules
	@variable(BUCR2model, bucrm_hydroGSpil[1:N_ZONES]>=0) # hydro energy schedules

	# Declaring decision variables for hourly dispatched and curtailed demand
	@variable(BUCR2model, bucrm_Demand[1:N_ZONES]>=0) # Hourly scheduled demand
	@variable(BUCR2model, bucrm_Demand_Curt[1:N_ZONES]>=0) # Hourly schedule demand
	@variable(BUCR2model, bucrm_OverGen[1:N_ZONES]>=0)

	# Declaring variables for transmission system
	@variable(BUCR2model, bucrm_voltAngle[1:N_ZONES]) #voltage angle at zone/bus n in t//
	@variable(BUCR2model, bucrm_powerFlow[1:N_ZONES, 1:M_ZONES]) #transmission Flow from zone n to zone m//

	# Defining the objective function that minimizes the overal cost of supplying electricity (=total variable, no-load, startup, and shutdown costs)
	@objective(BUCR2model, Min, sum(df_gens.IHRC_B1_HR[g]*fuelprice[g,day]*(bucrm_genOut_Block[g,1])
					   +df_gens.IHRC_B2_HR[g]*fuelprice[g,day]*(bucrm_genOut_Block[g,2])
					   +df_gens.IHRC_B3_HR[g]*fuelprice[g,day]*(bucrm_genOut_Block[g,3])
					   +df_gens.IHRC_B4_HR[g]*fuelprice[g,day]*(bucrm_genOut_Block[g,4])
					   +df_gens.IHRC_B5_HR[g]*fuelprice[g,day]*(bucrm_genOut_Block[g,5])
					   +df_gens.IHRC_B6_HR[g]*fuelprice[g,day]*(bucrm_genOut_Block[g,6])
					   +df_gens.IHRC_B7_HR[g]*fuelprice[g,day]*(bucrm_genOut_Block[g,7])
					   +(bucrm_TotGenVioP[g]*VIOLATION_PENALTY)
					   +(bucrm_TotGenVioN[g]*VIOLATION_PENALTY)
					   +(bucrm_MaxGenVioP[g]*VIOLATION_PENALTY)
					   +(bucrm_MinGenVioP[g]*VIOLATION_PENALTY)
					   +df_gens.NoLoadHR[g]*fuelprice[g,day]*(bucrm_genOnOff[g])
					   +((df_gens.HotStartU_FixedCost[g]
					   +(df_gens.HotStartU_HeatRate[g]*fuelprice[g,day]))*(bucrm_genStartUp[g]))
					   +df_gens.ShutdownCost[g]*(bucrm_genShutDown[g]) for g in 1:GENS)
					   +sum(df_peakers.IHRC_B1_HR[k]*fuelprice_peakers[k,day]*(bucrm_peakerOut_Block[k,1])
					   +df_peakers.IHRC_B2_HR[k]*fuelprice_peakers[k,day]*(bucrm_peakerOut_Block[k,2])
					   +df_peakers.IHRC_B3_HR[k]*fuelprice_peakers[k,day]*(bucrm_peakerOut_Block[k,3])
					   +df_peakers.IHRC_B4_HR[k]*fuelprice_peakers[k,day]*(bucrm_peakerOut_Block[k,4])
					   +df_peakers.IHRC_B5_HR[k]*fuelprice_peakers[k,day]*(bucrm_peakerOut_Block[k,5])
					   +df_peakers.IHRC_B6_HR[k]*fuelprice_peakers[k,day]*(bucrm_peakerOut_Block[k,6])
					   +df_peakers.IHRC_B7_HR[k]*fuelprice_peakers[k,day]*(bucrm_peakerOut_Block[k,7])
					   +df_peakers.NoLoadHR[k]*fuelprice_peakers[k,day]*(bucrm_peakerOnOff[k])
					   +((df_peakers.HotStartU_FixedCost[k]+(df_peakers.HotStartU_HeatRate[k]*fuelprice_peakers[k,day]))*(bucrm_peakerStartUp[k]))
					   +df_peakers.ShutdownCost[k]*(bucrm_peakerShutDown[k]) for k in 1:PEAKERS)
					   +sum((bucrm_Demand_Curt[n]*LOAD_SHED_PENALTY)+(bucrm_OverGen[n]*OVERGEN_PENALTY) for n=1:N_ZONES) )

	# Baseload Operation of nuclear units
	@constraint(BUCR2model, conNuckBaseLoad[g=1:GENS], bucrm_genOnOff[g]>=df_gens.Nuclear[g]) #
	@constraint(BUCR2model, conNuclearTotGenZone[n=1:N_ZONES], sum((bucrm_genOut[g]*map_gens[g,n]*df_gens.Nuclear[g]) for g=1:GENS) - bucr2_prepoc.nuclear[n] ==0)

	#Limits on generation of cogen units
	@constraint(BUCR2model, conCoGenBaseLoad[g=1:GENS], bucrm_genOnOff[g]>=df_gens.Cogen[g]) #
	@constraint(BUCR2model, conCoGenTotGenZone[n=1:N_ZONES], sum((bucrm_genOut[g]*map_gens[g,n]*df_gens.Cogen[g]) for g=1:GENS) - bucr2_prepoc.cogen[n] ==0)
	# Constraints representing technical limits of conventional generators
	#Status transition trajectory of
	@constraint(BUCR2model, conStartUpAndDn[g=1:GENS], (bucrm_genOnOff[g] - bucr2.gens.onoff_init[g] - bucrm_genStartUp[g] + bucrm_genShutDown[g])==0)
	# Max Power generation limit in Block 1
	@constraint(BUCR2model, conMaxPowBlock1[g=1:GENS],  bucrm_genOut_Block[g,1] <= df_gens.IHRC_B1_Q[g]*bucrm_genOnOff[g] )
	# Max Power generation limit in Block 2
	@constraint(BUCR2model, conMaxPowBlock2[g=1:GENS],  bucrm_genOut_Block[g,2] <= df_gens.IHRC_B2_Q[g]*bucrm_genOnOff[g] )
	# Max Power generation limit in Block 3
	@constraint(BUCR2model, conMaxPowBlock3[g=1:GENS],  bucrm_genOut_Block[g,3] <= df_gens.IHRC_B3_Q[g]*bucrm_genOnOff[g] )
	# Max Power generation limit in Block 4
	@constraint(BUCR2model, conMaxPowBlock4[g=1:GENS],  bucrm_genOut_Block[g,4] <= df_gens.IHRC_B4_Q[g]*bucrm_genOnOff[g] )
	# Max Power generation limit in Block 5
	@constraint(BUCR2model, conMaxPowBlock5[g=1:GENS],  bucrm_genOut_Block[g,5] <= df_gens.IHRC_B5_Q[g]*bucrm_genOnOff[g] )
	# Max Power generation limit in Block 6
	@constraint(BUCR2model, conMaxPowBlock6[g=1:GENS],  bucrm_genOut_Block[g,6] <= df_gens.IHRC_B6_Q[g]*bucrm_genOnOff[g] )
	# Max Power generation limit in Block 7
	@constraint(BUCR2model, conMaxPowBlock7[g=1:GENS],  bucrm_genOut_Block[g,7] <= df_gens.IHRC_B7_Q[g]*bucrm_genOnOff[g] )
	# Total Production of each generation equals the sum of generation from its all blocks
	@constraint(BUCR2model, conTotalGen[g=1:GENS],  sum(bucrm_genOut_Block[g,b] for b=1:BLOCKS) + bucrm_TotGenVioP[g]-bucrm_TotGenVioN[g]== bucrm_genOut[g])
	#Max power generation limit
	@constraint(BUCR2model, conMaxPow[g=1:GENS],  bucrm_genOut[g] - bucrm_MaxGenVioP[g] <= df_gens.MaxPowerOut[g]*bucrm_genOnOff[g])
	# Min power generation limit
	@constraint(BUCR2model, conMinPow[g=1:GENS],  bucrm_genOut[g] + bucrm_MinGenVioP[g]>= df_gens.MinPowerOut[g]*bucrm_genOnOff[g])
	#Up ramp rate limit
	@constraint(BUCR2model, conRampRateUp[g=1:GENS], (bucrm_genOut[g] - bucr2.gens.power_out_init[g] <=(df_gens.RampUpLimit[g]*bucr2.gens.onoff_init[g]) + (df_gens.RampStartUpLimit[g]*bucrm_genStartUp[g])))
	# Down ramp rate limit
	@constraint(BUCR2model, conRampRateDown[g=1:GENS], (bucr2.gens.power_out_init[g] - bucrm_genOut[g] <=(df_gens.RampDownLimit[g]*bucrm_genOnOff[g]) + (df_gens.RampShutDownLimit[g]*bucrm_genShutDown[g])))
	# Min Up Time limit with alternative formulation
	#The next twyo constraints enforce limits on binary commitment variables of slow and fast generators
	# scheduled slow units are forced to remain on, offline slow units remain off, and fast start units could change their commitment dependent on their MUT and MDT
	@constraint(BUCR2model, conCommitmentUB[g=1:GENS], (bucrm_genOnOff[g] <= bucr_commit_UB[g]))
	# if the generator is slow start and scheduled "on" in the SUCR,  is fixed by the following constraint
	@constraint(BUCR2model, conCommitmentLB[g=1:GENS], (bucrm_genOnOff[g] >= bucr_commit_LB[g]))

	# Constraints representing technical limits of peakers
	#Status transition trajectory of
	@constraint(BUCR2model, conStartUpAndDn_Peaker[k=1:PEAKERS], (bucrm_peakerOnOff[k] - bucr2.peakers.onoff_init[k] - bucrm_peakerStartUp[k] + bucrm_peakerShutDown[k])==0)
	# Max Power generation limit in Block 1
	@constraint(BUCR2model, conMaxPowBlock1_Peaker[k=1:PEAKERS],  bucrm_peakerOut_Block[k,1] <= df_peakers.IHRC_B1_Q[k]*bucrm_peakerOnOff[k] )
	# Max Power generation limit in Block 2
	@constraint(BUCR2model, conMaxPowBlock2_Peaker[k=1:PEAKERS],  bucrm_peakerOut_Block[k,2] <= df_peakers.IHRC_B2_Q[k]*bucrm_peakerOnOff[k] )
	# Max Power generation limit in Block 3
	@constraint(BUCR2model, conMaxPowBlock3_Peaker[k=1:PEAKERS],  bucrm_peakerOut_Block[k,3] <= df_peakers.IHRC_B3_Q[k]*bucrm_peakerOnOff[k] )
	# Max Power generation limit in Block 4
	@constraint(BUCR2model, conMaxPowBlock4_Peaker[k=1:PEAKERS],  bucrm_peakerOut_Block[k,4] <= df_peakers.IHRC_B4_Q[k]*bucrm_peakerOnOff[k] )
	# Max Power generation limit in Block 5
	@constraint(BUCR2model, conMaxPowBlock5_Peaker[k=1:PEAKERS],  bucrm_peakerOut_Block[k,5] <= df_peakers.IHRC_B5_Q[k]*bucrm_peakerOnOff[k] )
	# Max Power generation limit in Block 6
	@constraint(BUCR2model, conMaxPowBlock6_Peaker[k=1:PEAKERS],  bucrm_peakerOut_Block[k,6] <= df_peakers.IHRC_B6_Q[k]*bucrm_peakerOnOff[k] )
	# Max Power generation limit in Block 7
	@constraint(BUCR2model, conMaxPowBlock7_Peaker[k=1:PEAKERS],  bucrm_peakerOut_Block[k,7] <= df_peakers.IHRC_B7_Q[k]*bucrm_peakerOnOff[k] )
	# Total Production of each generation equals the sum of generation from its all blocks
	@constraint(BUCR2model, conTotalGen_Peaker[k=1:PEAKERS],  sum(bucrm_peakerOut_Block[k,b] for b=1:BLOCKS)>=bucrm_peakerOut[k])
	#Max power generation limit
	@constraint(BUCR2model, conMaxPow_Peaker[k=1:PEAKERS],  bucrm_peakerOut[k] <= df_peakers.MaxPowerOut[k]*bucrm_peakerOnOff[k])
	# Min power generation limit
	@constraint(BUCR2model, conMinPow_Peaker[k=1:PEAKERS],  bucrm_peakerOut[k] >= df_peakers.MinPowerOut[k]*bucrm_peakerOnOff[k])
	#Up ramp rate limit
	@constraint(BUCR2model, conRampRateUp_Peaker[k=1:PEAKERS], (bucrm_peakerOut[k] - bucr2.peakers.power_out_init[k] <=(df_peakers.RampUpLimit[k]*bucr2.peakers.onoff_init[k]) + (df_peakers.RampStartUpLimit[k]*bucrm_peakerStartUp[k])))
	# Down ramp rate limit
	@constraint(BUCR2model, conRampRateDown_Peaker[k=1:PEAKERS], (bucr2.peakers.power_out_init[k] - bucrm_peakerOut[k] <=(df_peakers.RampDownLimit[k]*bucrm_peakerOnOff[k]) + (df_peakers.RampShutDownLimit[k]*bucrm_peakerShutDown[k])))
	# Min Up Time limit with alternative formulation
	#The next twyo constraints enforce limits on binary commitment variables of slow and fast generators
	# scheduled slow units are forced to remain on, offline slow units remain off, and fast start units could change their commitment dependent on their MUT and MDT
	@constraint(BUCR2model, conCommitmentUB_Peaker[k=1:PEAKERS], (bucrm_peakerOnOff[k] <= bucr_commit_peaker_UB[k]))
	# if the generator is slow start and scheduled "on" in the SUCR,  is fixed by the following constraint
	@constraint(BUCR2model, conCommitmentLB_Peaker[k=1:PEAKERS], (bucrm_peakerOnOff[k] >= bucr_commit_peaker_LB[k]))

	# Renewable generation constraints
	@constraint(BUCR2model, conSolarLimit[n=1:N_ZONES], bucrm_solarG[n] + bucrm_solarGSpil[n] <= bucr2_prepoc.solar[n])
	@constraint(BUCR2model, conWindLimit[n=1:N_ZONES], bucrm_windG[n] + bucrm_windGSpil[n] <= bucr2_prepoc.wind[n])
	@constraint(BUCR2model, conHydroLimit[n=1:N_ZONES], bucrm_hydroG[n] + bucrm_hydroGSpil[n] <= bucr2_prepoc.hydro[n])

	# Constraints representing technical characteristics of storage units
	# the next three constraints fix the balancing charging/discharging/Idle status to their optimal outcomes as determined by SUCR
	@constraint(BUCR2model, conStorgChrgStatusFixed[p=1:STORG_UNITS], (bucrm_storgChrg[p]==sucr_to_bucr2.storg.chrg[p,hour]))
	@constraint(BUCR2model, conStorgDisChrgStatusFixed[p=1:STORG_UNITS], (bucrm_storgDisc[p]==sucr_to_bucr2.storg.disc[p,hour]))
	@constraint(BUCR2model, conStorgIdleStatusFixed[p=1:STORG_UNITS], (bucrm_storgIdle[p]==sucr_to_bucr2.storg.idle[p,hour]))

	# charging power limit
	@constraint(BUCR2model, conStrgChargPowerLimit[p=1:STORG_UNITS], (bucrm_storgChrgPwr[p] )<=df_storage.Power[p]*bucrm_storgChrg[p])
	# Discharging power limit
	@constraint(BUCR2model, conStrgDisChgPowerLimit[p=1:STORG_UNITS], (bucrm_storgDiscPwr[p])<=df_storage.Power[p]*bucrm_storgDisc[p])
	# State of charge at t
	@constraint(BUCR2model, conStorgSOC[p=1:STORG_UNITS], bucrm_storgSOC[p]==bucr2.storg.soc_init[p]-(bucrm_storgDiscPwr[p]/df_storage.TripEfficDown[p])+(bucrm_storgChrgPwr[p]*df_storage.TripEfficUp[p])-(bucrm_storgSOC[p]*df_storage.SelfDischarge[p]))
	# minimum energy limit
	@constraint(BUCR2model, conMinEnrgStorgLimi[p=1:STORG_UNITS], bucrm_storgSOC[p]>=0)
	# Maximum energy limit
	@constraint(BUCR2model, conMaxEnrgStorgLimi[p=1:STORG_UNITS], bucrm_storgSOC[p]<=(df_storage.Power[p]/df_storage.PowerToEnergRatio[p]))

	# Constraints representing transmission grid capacity constraints
	# DC Power Flow Calculation
	#@constraint(BUCR2model, conDCPowerFlowPos[n=1:N_ZONES, m=1:N_ZONES], bucrm_powerFlow[n,m]-(TranS[n,m]*(bucrm_voltAngle[n]-bucrm_voltAngle[m])) ==0)
	@constraint(BUCR2model, conDCPowerFlowNeg[n=1:N_ZONES, m=1:N_ZONES], bucrm_powerFlow[n,m]+bucrm_powerFlow[m,n]==0)
	# Tranmission flow bounds (from n to m and from m to n)
	@constraint(BUCR2model, conPosFlowLimit[n=1:N_ZONES, m=1:N_ZONES], bucrm_powerFlow[n,m]<=TranC[n,m])
	@constraint(BUCR2model, conNegFlowLimit[n=1:N_ZONES, m=1:N_ZONES], bucrm_powerFlow[n,m]>=-TranC[n,m])
	# Voltage Angle bounds and reference point
	#@constraint(BUCR2model, conVoltAnglUB[n=1:N_ZONES], bucrm_voltAngle[n]<=π)
	#@constraint(BUCR2model, conVoltAnglLB[n=1:N_ZONES], bucrm_voltAngle[n]>=-π)
	#@constraint(BUCR2model, conVoltAngRef, bucrm_voltAngle[1]==0)

	# Demand-side Constraints
	@constraint(BUCR2model, conDemandLimit[n=1:N_ZONES], bucrm_Demand[n]+ bucrm_Demand_Curt[n] == bucr2_prepoc.demand[n])

	# Demand Curtailment and wind generation limits
	@constraint(BUCR2model, conDemandCurtLimit[n=1:N_ZONES], bucrm_Demand_Curt[n] <= LOAD_SHED_MAX);
	@constraint(BUCR2model, conOverGenLimit[n=1:N_ZONES], bucrm_OverGen[n] <= OVERGEN_MAX);

	# System-wide Constraints
	#nodal balance constraint
	@constraint(BUCR2model, conNodBalanc[n=1:N_ZONES], sum((bucrm_genOut[g]*map_gens[g,n]) for g=1:GENS) + sum((bucrm_peakerOut[k]*map_peakers[k,n]) for k=1:PEAKERS) + sum((bucrm_storgDiscPwr[p]*map_storage[p,n]) for p=1:STORG_UNITS) - sum((bucrm_storgChrgPwr[p]*map_storage[p,n]) for p=1:STORG_UNITS) +bucrm_solarG[n] +bucrm_windG[n] +bucrm_hydroG[n] - bucrm_OverGen[n]- bucrm_Demand[n] == sum(bucrm_powerFlow[n,m] for m=1:M_ZONES))
	# Minimum up reserve requirement
	#    @constraint(BUCR2model, conMinUpReserveReq[t=1:N_Hrs_BUCR], sum(genResUp[g,t] for g=1:GENS) + sum(storgResUp[p,t] for p=1:STORG_UNITS) >= Reserve_Req_Up[t] )

	# Minimum down reserve requirement
	#    @constraint(BUCR2model, conMinDnReserveReq[t=1:N_Hrs_BUCR], sum(genResDn[g,t] for g=1:GENS) + sum(storgResDn[p,t] for p=1:STORG_UNITS) >= Reserve_Req_Dn[t] )

	t2_BUCR2model = time_ns()
	time_BUCR2model = (t2_BUCR2model -t1_BUCR2model)/1.0e9;
	@info "BUCR2model for day: $day and hour $(hour+INIT_HR_SUCR) setup executed in (s): $time_BUCR2model";

	open(".//outputs//time_performance.csv", FILE_ACCESS_APPEND) do io
			writedlm(io, hcat("BUCR2model", time_BUCR2model, "day: $day",
					"hour $(hour+INIT_HR_SUCR)", "", "Model Setup"), ',')
	end; # closes file

	# solve the First WAUC model (BUCR)
	JuMP.optimize!(BUCR2model)
	open(".//outputs//objective_values_v76.csv", FILE_ACCESS_APPEND) do io
			writedlm(io, hcat("BUCR2model", "day: $day",
					"hour $(hour+INIT_HR_SUCR)", "", JuMP.objective_value(BUCR2model)), ',')
	end;

	#compute_conflict!(BUCR2model)
	#if MOI.get(BUCR2model, MOI.ConflictStatus()) != MOI.CONFLICT_FOUND
	#    error("No conflict could be found for an infeasible model.")
	#end

	# Both constraints should participate in the conflict.
	#MOI.get(BUCR2model, MOI.ConstraintConflictStatus(), conTotalGen)
	#MOI.get(BUCR2model, MOI.ConstraintConflictStatus(), c2)

	# Pricing general results in the terminal window
	println("Objective value: ", JuMP.objective_value(BUCR2model))

	println("------------------------------------")
	println("------- BAUC2 OBJECTIVE VALUE -------")
	println("Objective value for day ", day, " and hour ", hour+INIT_HR_SUCR,
			" is: ", JuMP.objective_value(BUCR2model))
	println("------------------------------------")
	println("------- BAUC2 PRIMAL STATUS -------")
	println(primal_status(BUCR2model))
	println("------------------------------------")
	println("------- BAUC2 DUAL STATUS -------")
	println(JuMP.dual_status(BUCR2model))
	println("Day: ", day, " and hour ", hour+INIT_HR_SUCR, " solved")
	println("---------------------------")
	println("BUCR2model Number of variables: ", JuMP.num_variables(BUCR2model))
	@info "BUCR2model Number of variables: " JuMP.num_variables(BUCR2model)

	open(".//outputs//time_performance.csv", FILE_ACCESS_APPEND) do io
			writedlm(io, hcat("BUCR2model", JuMP.num_variables(BUCR2model), "day: $day",
					"hour $(hour+INIT_HR_SUCR)", "", "Variables"), ',')
	end;

	@debug "BUCR2model for day: $day and hour $(hour+INIT_HR_SUCR) optimized executed in (s): $(solve_time(BUCR2model))";

	open(".//outputs//time_performance.csv", FILE_ACCESS_APPEND) do io
			writedlm(io, hcat("BUCR2model", solve_time(BUCR2model), "day: $day",
					"hour $(hour+INIT_HR_SUCR)", "", "Model Optimization"), ',')
	end; # closes file

	## Write the optimal outcomes into spreadsheets
	# TODO: Later we need to include a variable for day so the cell number in
	# which the results are printed is updated accordingly

	# Write the conventional generators' schedules
	t1_write_BUCR2model_results = time_ns()
	open(".//outputs//BUCR_GenOutputs.csv", FILE_ACCESS_APPEND) do io
		for g=1:GENS
			writedlm(io, hcat(day, hour+INIT_HR_SUCR, g, df_gens.UNIT_ID[g],
				df_gens.MinPowerOut[g], df_gens.MaxPowerOut[g],
				JuMP.value.(bucrm_genOut[g]), JuMP.value.(bucrm_genOnOff[g]),
				JuMP.value.(bucrm_genShutDown[g]),
				JuMP.value.(bucrm_genStartUp[g]),
				JuMP.value.(bucrm_TotGenVioP[g]),
				JuMP.value.(bucrm_TotGenVioN[g]),
				JuMP.value.(bucrm_MaxGenVioP[g]),
				JuMP.value.(bucrm_MinGenVioP[g])), ',')
		end # ends the loop
	end; # closes file

	# Write the peakers' schedules
	open(".//outputs//BUCR_PeakerOutputs.csv", FILE_ACCESS_APPEND) do io
		for k=1:PEAKERS
			writedlm(io, hcat(day, hour+INIT_HR_SUCR, k, df_peakers.UNIT_ID[k],
				df_peakers.MinPowerOut[k], df_peakers.MaxPowerOut[k],
				JuMP.value.(bucrm_peakerOut[k]),
				JuMP.value.(bucrm_peakerOnOff[k]),
				JuMP.value.(bucrm_peakerShutDown[k]),
				JuMP.value.(bucrm_peakerStartUp[k])), ',')
		end # ends the loop
	end; # closes file

	# Writing storage units' optimal schedules into CSV file
	open(".//outputs//BUCR_StorageOutputs.csv", FILE_ACCESS_APPEND) do io
		for p=1:STORG_UNITS
			writedlm(io, hcat(day, hour+INIT_HR_SUCR, p, df_storage.Name[p],
					df_storage.Power[p],
					df_storage.Power[p]/df_storage.PowerToEnergRatio[p],
					JuMP.value.(bucrm_storgChrg[p]),
					JuMP.value.(bucrm_storgDisc[p]),
					JuMP.value.(bucrm_storgIdle[p]),
					JuMP.value.(bucrm_storgChrgPwr[p]),
					JuMP.value.(bucrm_storgDiscPwr[p]),
					JuMP.value.(bucrm_storgSOC[p]) ), ',')
		end # ends the loop
	end; # closes file

	# Writing the transmission flow schedules in CSV file
	open(".//outputs//BUCR_TranFlowOutputs.csv", FILE_ACCESS_APPEND) do io
		for n=1:N_ZONES, m=1:M_ZONES
			writedlm(io, hcat(day, hour+INIT_HR_SUCR, n, m,
				JuMP.value.(bucrm_powerFlow[n,m]), TranC[n,m]), ',')
		end # ends the loop
	end; # closes file

	# Writing the curtilment, overgeneration, and spillage outcomes in CSV file
	open(".//outputs//BUCR_Curtail.csv", FILE_ACCESS_APPEND) do io
		for n=1:N_ZONES
		   writedlm(io, hcat(day, hour+INIT_HR_SUCR, n,
			   JuMP.value.(bucrm_OverGen[n]), JuMP.value.(bucrm_Demand_Curt[n]),
			   JuMP.value.(bucrm_windGSpil[n]), JuMP.value.(bucrm_solarGSpil[n]),
			   JuMP.value.(bucrm_hydroGSpil[n])), ',')
		end # ends the loop
	end; # closes file
	t2_write_BUCR2model_results = time_ns()

	time_write_BUCR2model_results = (t2_write_BUCR2model_results -t1_write_BUCR2model_results)/1.0e9;
	@info "Write BUCR2model results for day $day and hour $(hour+INIT_HR_SUCR) executed in (s): $time_write_BUCR2model_results";

	open(".//outputs//time_performance.csv", FILE_ACCESS_APPEND) do io
			writedlm(io, hcat("BUCR2model", time_write_BUCR2model_results, "day: $day",
					"hour $(hour+INIT_HR_SUCR)", "", "Write CSV files"), ',')
	end; #closes file

	# Initilization of the next UC Run
	# Setting initial values for conventional generators in BUCR2 (next hour),
	# FUCR, and BUCR1
	for g=1:GENS
		# set the initiali values to be fed to the next hour BUCR2
		bucr2.gens.onoff_init[g] = JuMP.value.(bucrm_genOnOff[g]); #
		bucr2.gens.power_out_init[g] = JuMP.value.(bucrm_genOut[g]);
		# Set the initial values fed to the next FUCR and BUCR1

		if hour==(24-INIT_HR_SUCR+INIT_HR_FUCR)
			bucr1.gens.onoff_init[g] = JuMP.value.(bucrm_genOnOff[g]);
			bucr1.gens.power_out_init[g] = JuMP.value.(bucrm_genOut[g]);
			fucr.gens.onoff_init[g] = JuMP.value.(bucrm_genOnOff[g]);
			fucr.gens.power_out_init[g] = JuMP.value.(bucrm_genOut[g]);
		end
	end

	# Setting initial values for peakers in BUCR2 (next hour), FUCR, and BUCR1
	for g=1:PEAKERS
		# set the initiali values to be fed to the next hour BUCR2
		bucr2.peakers.onoff_init[g] = round(JuMP.value.(bucrm_peakerOnOff[g]));
		bucr2.peakers.power_out_init[g] = JuMP.value.(bucrm_peakerOut[g]);
		# Set the initial values fed to the next FUCR and BUCR1

		if hour==(24-INIT_HR_SUCR+INIT_HR_FUCR)
			bucr1.peakers.onoff_init[g] = round(JuMP.value.(bucrm_peakerOnOff[g]));
			bucr1.peakers.power_out_init[g] = JuMP.value.(bucrm_peakerOut[g]);
			fucr.peakers.onoff_init[g] = round(JuMP.value.(bucrm_peakerOnOff[g]));
			fucr.peakers.power_out_init[g] = JuMP.value.(bucrm_peakerOut[g]);
		end
	end

	for p=1:STORG_UNITS
		# Set the initial values to be fed to the next hour BUCR2
		bucr2.storg.soc_init[p]=JuMP.value.(bucrm_storgSOC[p]);
		# set the initiali values to be fed to the next hour BUCR2
		if hour==(24-INIT_HR_SUCR+INIT_HR_FUCR)
			bucr1.storg.soc_init[p] = JuMP.value.(bucrm_storgSOC[p]);
			fucr.storg.soc_init[p] = JuMP.value.(bucrm_storgSOC[p]);
		end
	end

	#Update the up and down times for individual  generators
	for g=1:GENS
		# Update the total up-time or down-time for each individual conventional
		# generator to be fed into BUCR2
		if JuMP.value.(bucrm_genStartUp[g]) > 0.5
			bucr2.gens.uptime_init[g]= 1;
			bucr2.gens.dntime_init[g] = 0;
		elseif JuMP.value.(bucrm_genShutDown[g]) > 0.5
			bucr2.gens.uptime_init[g]= 0;
			bucr2.gens.dntime_init[g]= 1;
		else
			if JuMP.value.(bucrm_genOnOff[g]) > 0.5
				bucr2.gens.uptime_init[g]=bucr2.gens.uptime_init[g]+1;
				bucr2.gens.dntime_init[g]= 0;
			else
				bucr2.gens.uptime_init[g]= 0;
				bucr2.gens.dntime_init[g]= bucr2.gens.dntime_init[g]+1;
			end
		end
		#Update the total up and down times to be fed into FUCR and BUCR1
		if hour == (INIT_HR_SUCR-INIT_HR_FUCR)
			bucr1.gens.uptime_init[g]= bucr2.gens.uptime_init[g];
			bucr1.gens.dntime_init[g]= bucr2.gens.dntime_init[g];
			fucr.gens.uptime_init[g]= bucr2.gens.uptime_init[g];
			fucr.gens.dntime_init[g]= bucr2.gens.dntime_init[g];
		end
	end

	for g = 1:PEAKERS
		# Update the total up-time or down-time for each individual
		# peakers to be fed into BUCR2
		if (JuMP.value.(bucrm_peakerStartUp[g])) > 0
			bucr2.peakers.uptime_init[g]= 1;
			bucr2.peakers.dntime_init[g] = 0;
		elseif (JuMP.value.(bucrm_peakerShutDown[g])) > 0
			bucr2.peakers.uptime_init[g]= 0;
			bucr2.peakers.dntime_init[g]= 1;
		else
			if (JuMP.value.(bucrm_peakerOnOff[g])) > 0
				bucr2.peakers.uptime_init[g]= bucr2.peakers.uptime_init[g]+1;
				bucr2.peakers.dntime_init[g]= 0;
			else
				bucr2.peakers.uptime_init[g]= 0;
				bucr2.peakers.dntime_init[g]= bucr2.peakers.dntime_init[g]+1;
			end
		end
		#Update the total up and down times to be fed into FUCR and BUCR1
		if hour == (INIT_HR_SUCR-INIT_HR_FUCR)
			bucr1.peakers.uptime_init[g]= bucr2.peakers.uptime_init[g];
			bucr1.peakers.dntime_init[g]= bucr2.peakers.dntime_init[g];
			fucr.peakers.uptime_init[g]= bucr2.peakers.uptime_init[g];
			fucr.peakers.dntime_init[g]= bucr2.peakers.dntime_init[g];
		end
	end

	return fucr, bucr1, bucr2;

end; # end function fucr2_model

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
"""
    MVCalc(μstar,μ,Σ)

Calculate the std and weights of a portfolio (with mean return μstar) on MVF of risky assets.

# Remark
- Only (λ,δ) and thus (w,stdRp) depend on μstar. We could therefore speed up the computations a bit
by doing the loop over different μstar values inside the function (and thus not recalculate Σ_1,a,b,c).
"""
function MVCalc(μstar,μ,Σ)
    n    = length(μ)
    Σ_1  = inv(Σ)
    a    = μ'Σ_1*μ
    b    = μ'Σ_1*ones(n)
    c    = ones(n)'Σ_1*ones(n)
    λ    = (c*μstar - b)/(a*c-b^2)
    δ    = (a-b*μstar)/(a*c-b^2)
    w    = Σ_1 *(μ*λ.+δ)
    StdRp = sqrt(w'Σ*w)
    return StdRp,w
end

"""
    MVCalcRf(μstar,μ,Σ,Rf)

Calculate the std and portfolio weights of a portfolio (with a given mean, μstar) on MVF of (risky assets,riskfree)
"""
function MVCalcRf(μstar,μ,Σ,Rf)
    μe    = μ .- Rf
    Σ_1   = inv(Σ)
    w     = (μstar-Rf)/(μe'Σ_1*μe) * Σ_1*μe
    StdRp = sqrt(w'Σ*w)
    return StdRp,w                    #std and portfolio weights
end

"""
    MVTangencyP(μ,Σ,Rf)

Calculate the tangency portfolio
"""
function MVTangencyP(μ,Σ,Rf)           #calculates the tangency portfolio
    n    = length(μ)
    μe   = μ .- Rf                    #expected excess returns
    Σ_1  = inv(Σ)
    w    = Σ_1 *μe/(ones(n)'Σ_1*μe)
    muT  = w'μ + (1-sum(w))*Rf
    StdT = sqrt(w'Σ*w)
    return w,muT,StdT                  #portolio weights, mean and std
end

ext: jl

max_stars_repo_name: snowdj/FinancialTheoryMSc

max_issues_repo_name: snowdj/FinancialTheoryMSc

max_forks_repo_name: snowdj/FinancialTheoryMSc


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SCS_GPU_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SCS_GPU")
JLLWrappers.@generate_main_file("SCS_GPU", UUID("af6e375f-46ec-5fa0-b791-491b0dfa44a4"))
end  # module SCS_GPU_jll


text:
using Pkg
const LONGER_TESTS = false
const GROUP = get(ENV, "GROUP", "All")

using OperatorLearning
using Test
using Random

Random.seed!(0)

if GROUP == "All" || GROUP == "Core"
@testset "FourierLayer" begin
    include("fourierlayer.jl")
end

@testset "DeepONet" begin
    include("deeponet.jl")
end

@testset "Weights" begin
    include("complexweights.jl")
end
end

if GROUP == "GPU"
    # Add GPU Tests Here
end


text:
"""
	2D IMPLEMENTATION
"""

# crosscorr( f, g ) = iFFT( conj( FFT( f ) ) .* FFT( g ) )

function crossCorrelation!( ::FFT, cmat::A{T,2}, padf::A{C{T},2}, padg::A{C{T},2}, plan, iplan
                          ) where {T<:AbstractFloat}

    plan * padf;
	plan * padg;
    @inbounds @simd for e in 1:length(padf)
        padf[e] = conj( padf[e] ) * padg[e];
    end

	iplan * padf;
	@inbounds @simd for e in 1:length(padf)
		cmat[e] = real( padf[e] )
	end
end

function PIV_2D( ::FFT, img1::A{<:Real,2}, img2::A{<:Real,2},
                        IA::III, SM::III, overlap::III, mpass::I, width::I,
						peak::S, sigNoise::S, filtFun::Function, threshold::F;
 					    corrType=Float32, vfType=Float32 )

    step   = IA  .-  overlap;
	imsize = ( size(img1,1), size(img1,2), size(img1,3) )
    vfSize = length.( StepRange.( IA, step, imsize ) );

    U  = zeros( vfType, vfSize[1:2] );
    V  = zeros( vfType, vfSize[1:2] );
    SN = zeros( vfType, vfSize[1:2] );

	ignoreSN = sigNoise == "";
	SM = ( SM[1], SM[2], 0 );

    for mp in mpass:-1:1

        last_mp = ( mp == 1 );

        # Scaling IA, SM and step (overlap) to multi-pass.
        IA_mp    =  IA  .* mp;
        SM_mp    =  SM  .* mp;
        step_mp  = step .* mp;
		IAranges = StepRange.( IA_mp, step_mp, imsize )

        # Initializing variables for FFT cross-Correlation.
        csize   = 2 .* ( IA_mp .+ SM_mp )
        shifts  = div.( csize, 2 ) .+ SM_mp .- 1;
        cmatrix = zeros( corrType, csize[1:2] );
        shifted = zeros( corrType, csize[1:2] );
        pads    = zeros( Complex{corrType}, csize[1:2] );
        padi    = zeros( Complex{corrType}, csize[1:2] );
        plan    = FFTW.plan_fft!(  pads );
        iplan   = FFTW.plan_ifft!( pads );

		# Nested loop establishing IA and UV coordinates
		vfx = [ 1-mp, 0 ]
		for x2 in IAranges[2];     x1 = x2-IA_mp[2]+1; vfx .+= mp; vfy = [ 1-mp, 0 ];
			for y2 in IAranges[1]; y1 = y2-IA_mp[1]+1; vfy .+= mp;

				# 0-. Filtering undesirable IAs, ex background
				interr = view( img1, y1:y2, x1:x2 )

				if filtFun( interr ) < threshold
					SN[ vfy[1]:vfy[2], vfx[1]:vfx[2] ] .= -1.0
					continue;
				end

			    # 1.1-. Previously computed displacements shift the search area during multipass
			    u = round(Int32, U[ vfy[1], vfx[1] ]);
			    v = round(Int32, V[ vfy[1], vfx[1] ]);

				# 1.2-. Search volume coordinates after shifting.
				scoords, offset = SAcoords( (y1,x1,1), (y2,x2,1), (u,v,0), SM_mp, imsize );

				# 1.3-. Copying IA and SA into the padded arrays for FFT
				setTo0!( pads, padi )
			    putWithinPadded!( pads, img2, 0.0, offset, scoords );
			    putWithinPadded!( padi, img1, 0.0, (0,0,0), (y1,y2,x1,x2,1,1) );

           		# 2-. cross-correlation
		        crossCorrelation!( FFT(), cmatrix, pads, padi, plan, iplan )
		        Base.circshift!( shifted, cmatrix, shifts[1:2] )

			    # 3.1-. Extracting the displacement
			    ( r, c ) = approxTranslation( shifted, peak, Val(last_mp) )

				# 3.2 -. Computing signal to noise ration, if needed
	 			if ( !ignoreSN && last_mp )
			    	SN[ vfy[1]:vfy[2], vfx[1]:vfx[2] ] .= SNRatio( shifted, sigNoise, width )
				end

			    # 4-. Updating U, V matrices
			    U[ vfy[1]:vfy[2], vfx[1]:vfx[2] ] .= u - r;
			    V[ vfy[1]:vfy[2], vfx[1]:vfx[2] ] .= v - c;
		end end
    end # Multi-pass loop

    return U, V, SN
end

"""
    3D IMPLEMENTATION
"""

# crosscorr( f, g ) = iFFT( conj( FFT( f ) ) .* FFT( g ) )
function crossCorrelation!( ::FFT, cmat::A{T,3}, padf::A{C{T},3}, padg::A{C{T},3}, plan, iplan
                          ) where {T<:AbstractFloat}
    plan * padf;
    plan * padg;
    @inbounds @simd for e in 1:length(padf)
        padf[e] = conj( padf[e] ) * padg[e];
    end

    iplan * padf;
    @inbounds @simd for e in 1:length(padf)
        cmat[e] = real( padf[e] )
    end
end

function PIV_3D( ::FFT, vol1::A{T,3}, vol2::A{T,3},
                        IA::III, SM::III, overlap::III, mpass::I, width::I,
                        peak::S, sigNoise::S, filtFun::Function, threshold::F ;
                        corrType=Float32, vfType=Float32
               ) where {T<:Real}

	# Calculating vector field size.
    step   = IA  .-  overlap;
    vfSize = length.( StepRange.( IA, step, size(vol1) ) );

    U  = zeros( vfType, vfSize );
    V  = zeros( vfType, vfSize );
    W  = zeros( vfType, vfSize );
    SN = zeros( vfType, vfSize );

	ignoreSN = sigNoise == "";

    for mp in mpass:-1:1

        last_mp = ( mp == 1 );

        # Scaling IA, SM and step (overlap) to multi-pass iteration.
        IA_mp    =  IA  .* mp;
        SM_mp    =  SM  .* mp;
        step_mp  = step .* mp;
		IAranges = StepRange.( IA_mp, step_mp, size(vol1) )

        # Initialize cross-correlation variables once per multi-pass iteration
        csize   = 2 .* ( IA_mp .+ SM_mp )
        shifts  = div.( csize, 2 ) .+ SM_mp .- 1;
        cmatrix = zeros( corrType, csize );
        shifted = zeros( corrType, csize );
        padi    = zeros( Complex{corrType}, csize );
        pads    = zeros( Complex{corrType}, csize );
        plan    = FFTW.plan_fft!(  pads );
        iplan   = FFTW.plan_ifft!( pads );

		# Nested loop establishing IV and UVW coordinates
		vfz = [ 1-mp, 0 ]
		for z2 in IAranges[3];         z1 = z2-IA_mp[3]+1; vfz .+= mp; vfx = [ 1-mp, 0 ];
			for x2 in IAranges[2];     x1 = x2-IA_mp[2]+1; vfx .+= mp; vfy = [ 1-mp, 0 ];
				for y2 in IAranges[1]; y1 = y2-IA_mp[1]+1; vfy .+= mp;

					# 0-. Filtering undesirable IAs, ex background
					interr = view( vol1, y1:y2, x1:x2, z1:z2 )

					if filtFun( interr ) < threshold
						SN[ vfy[1]:vfy[2], vfx[1]:vfx[2], vfz[1]:vfz[2] ] .= NaN32
						continue;
					end

				    # 1-. Previously computed displacements shift the search area during multipass
				    u = round(Int32, U[ vfy[1], vfx[1], vfz[1] ]);
				    v = round(Int32, V[ vfy[1], vfx[1], vfz[1] ]);
				    w = round(Int32, W[ vfy[1], vfx[1], vfz[1] ]);

				    # 1.2-. Search volume coordinates after shifting.
					scoords, off = SAcoords( (y1,x1,z1), (y2,x2,z2), (u,v,w), SM_mp, size(vol1) );

					# 1.3-. Copying IV and SV into the padded arrays for FFT
					setTo0!( pads )
					setTo0!( padi )
				    putWithinPadded!( pads, vol2, 0.0, off, scoords );
				    putWithinPadded!( padi, vol1, 0.0, (0,0,0), (y1,y2,x1,x2,z1,z2) );

					# 2-. Cross-correlation
					crossCorrelation!( FFT(), cmatrix, pads, padi, plan, iplan )
		            Base.circshift!( shifted, cmatrix, shifts );

				    # 3-. Extracting the displacement
				    ( r, c, z ) = approxTranslation( shifted, peak, Val(last_mp) )

					# 3.2 -. Computing signal to noise ration, if needed
		 			if ( !ignoreSN && last_mp )
						snres = SNRatio( shifted, sigNoise, width )
				    	SN[ vfy[1]:vfy[2], vfx[1]:vfx[2], vfz[1]:vfz[2] ] .= snres
					end

				    # 4-. Updating U, V, W matrices
				    U[ vfy[1]:vfy[2], vfx[1]:vfx[2], vfz[1]:vfz[2] ] .= u - r;
				    V[ vfy[1]:vfy[2], vfx[1]:vfx[2], vfz[1]:vfz[2] ] .= v - c;
				    W[ vfy[1]:vfy[2], vfx[1]:vfx[2], vfz[1]:vfz[2] ] .= w - z;
		end	end	end
    end # Multi-pass loop

    return U, V, W, SN
end

""" FFT out of place CROSS CORRELATION """

function crossCorrelation( ::FFT, g::A{<:Real,N}, f::A{<:Real,N}; corrType=Float32, shift=false ) where {N}

    csize = size( f ) .+ size( g );

    corr  = zeros( corrType, csize );
    padf  = zeros( Complex{corrType}, csize );
    padg  = zeros( Complex{corrType}, csize );
    plan  =  plan_fft!( padf );
    iplan = plan_ifft!( padf );

    putWithinPadded!( padf, f, 0.0, (0,0,0), (1,size(f,1),1,size(f,2),1,size(f,3)) );
    putWithinPadded!( padg, g, 0.0, (0,0,0), (1,size(g,1),1,size(g,2),1,size(f,3)) );

	crossCorrelation!( FFT(), corr, padf, padg, plan, iplan )

    if shift
        shifted = copy( corr );
        shifts  = div.( csize, 2 ) .+ div.( size(f) .- size(g), 2 ) .- 1;
        Base.circshift!( corr, shifted, shifts );
    end

   return corr
end

max_issues_repo_name: ErichZimmer/quickPIV

max_forks_repo_name: ErichZimmer/quickPIV


text:
using Images

function imwarp(
    pts::Array{Float64,2},
    img::Image,
    nCols::Int64,
    nRows::Int64,
  )

  newImg::Image = copyproperties(img, Array(typeof(img[1,1]), nCols, nRows))
  x::Float64 = 0.0
  y::Float64 = 0.0
  x0::Int64 = 0
  y0::Int64 = 0
  x1::Int64 = 0
  y1::Int64 = 0
  xx::Float64 = 0.0
  yy::Float64 = 0.0
  for i::Int64=1:nCols
    for j::Int64=1:nRows
      x = pts[i+(j-1)*nRows,1]
      y = pts[i+(j-1)*nRows,2]

      x0 = floor(Int64, x)
      y0 = floor(Int64, y)
      x1 = ceil(Int64, x)
      y1 = ceil(Int64, y)

      if (x0 < 0 || y0 < 0 || x1 >= nCols || y1 >= nRows)
        newImg.data[i,j] = 0*img[1,1]
        continue
      end

      xx = x-x0
      yy = y-y0
      newImg.data[i, j] =
        img.data[x0+1, y0+1]*(1-xx)*(1-yy)+
        img.data[x0+1, y1+1]*xx*(1-yy)+
        img.data[x1+1, y0+1]*(1-xx)*yy+
        img.data[x1+1, y1+1]*xx*yy
    end
  end
  newImg
end


text:
import Clang
import Cbc_jll
import CoinUtils_jll

FILES = [
    (CoinUtils_jll.artifact_dir, "Coin_C_defines.h"),
    (Cbc_jll.artifact_dir, "Cbc_C_Interface.h"),
]

for (dir, file) in FILES
    cp(
        joinpath(dir, "include", "coin", file),
        joinpath(@__DIR__, file);
        force = true,
    )
end

GEN_DIR = joinpath(@__DIR__, "..", "src", "gen")

wc = Clang.init(
    headers = map(i -> joinpath(@__DIR__, i[2]), FILES),
    output_file = joinpath(GEN_DIR, "libcbc_api.jl"),
    common_file = joinpath(GEN_DIR, "libcbc_common.jl"),
    header_wrapped = (root, current) -> root == current,
    header_library = x -> "libcbcsolver",
    clang_diagnostics = true,
)

run(wc)

rm(joinpath(GEN_DIR, "LibTemplate.jl"))


text:
export reconstruction_direct_2d, reconstruction_direct_3d

"""
    reconstruction_direct(acqData::AcquisitionData, reconSize::NTuple{D,Int64}, weights::Vector{Vector{Complex{<:AbstractFloat}}}, correctionMap::Array{Complex{<:AbstractFloat}}=Complex{<:AbstractFloat}[])

Performs a direct Fourier-based image reconstruction of AcquisitionData

input:
  `acqData::AcquisitionData`            - AcquisitionData object
  `reconSize::NTuple{D,Int64}`              - size of image to reconstruct
  `weights::Vector{Vector{Complex{<:AbstractFloat}}}` - sampling density of the trajectories in acqData
  (`correctionMap::Array{Complex{<:AbstractFloat}}`)  - fieldmap for the correction of off-resonance effects

"""
function reconstruction_direct(acqData::AcquisitionData{T}
                                  , reconSize::NTuple{D,Int64}
                                  , weights::Vector{Vector{Complex{T}}}
                                  , correctionMap::Array{Complex{T}}=Complex{T}[]) where {D,T}

  encDims = ndims(trajectory(acqData))
  if encDims!=D
    error("reco-dimensionality $D and encoding-dimensionality $(encDims) do not match")
  end

  numContr, numChan, numSl = numContrasts(acqData), numChannels(acqData), numSlices(acqData)
  Ireco = zeros(Complex{T}, prod(reconSize), numSl, numContr, numChan)

  p = Progress(numSl*numChan*numContr, 1, "Direct Reconstruction...")

  for i = 1:numSl
    F = encodingOps_simple(acqData, reconSize, slice=i, correctionMap=correctionMap)
    for k = 1:numContr
      for j = 1:numChan
        kdata = kData(acqData,k,j,i) .* (weights[k].^2)
        I = adjoint(F[k]) * kdata

        if isCircular( trajectory(acqData, k) )
          circularShutter!(reshape(I, reconSize), 1.0)
        end
        Ireco[:,i,k,j] = I

        next!(p)
      end
    end
  end

  Ireco = reshape(Ireco, volumeSize(reconSize, numSl)..., numContr, numChan)
  # if D==2
  #   # 2d reconstruction
  #   Ireco = reshape(Ireco, reconSize[1], reconSize[2], numSl, numContr, numChan)
  # else
  #   # 3d reconstruction
  #   Ireco = reshape(Ireco, reconSize[1], reconSize[2], reconSize[3], numContr, numChan)
  # end

  return makeAxisArray(Ireco, acqData)
end


text:
using JuMP
using LinearAlgebra

function constraint_def_abs_flow_values(pm::AbstractPowerModel, l::Int; nw::Int=pm.cnw)
    branch = ref(pm, nw, :branch, l)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (l, f_bus, t_bus)
    t_idx = (l, t_bus, f_bus)

    cref_pf,cref_pt=nothing,nothing
    if haskey( var(pm,nw), :p )
      p_fr = var(pm, nw,:p, f_idx)
      p_to = var(pm, nw,:p, t_idx)

      pd_fm = var(pm,nw, :pd_br)[f_idx,0]
      pd_fp = var(pm,nw, :pd_br)[f_idx,1]
      pd_tm = var(pm,nw, :pd_br)[t_idx,0]
      pd_tp = var(pm,nw, :pd_br)[t_idx,1]

      pt_fm = var(pm,nw, :pt_br)[f_idx,0]
      pt_fp = var(pm,nw, :pt_br)[f_idx,1]
      pt_tm = var(pm,nw, :pt_br)[t_idx,0]
      pt_tp = var(pm,nw, :pt_br)[t_idx,1]

      cref_pf=JuMP.@constraint(pm.model, p_fr - pd_fm + pd_fp + pt_fm - pt_fp == 0)
      cref_pt=JuMP.@constraint(pm.model, p_to - pd_tm + pd_tp + pt_tm - pt_tp == 0)
    end

    cref_qf,cref_qt=nothing,nothing
    if haskey( var(pm,nw), :q )
      q_fr = var(pm, nw, :q, f_idx)
      q_to = var(pm, nw, :q, t_idx)

      qd_fm = var(pm,nw, :qd_br)[f_idx,0]
      qd_fp = var(pm,nw, :qd_br)[f_idx,1]
      qd_tm = var(pm,nw, :qd_br)[t_idx,0]
      qd_tp = var(pm,nw, :qd_br)[t_idx,1]

      qt_fm = var(pm,nw, :qt_br)[f_idx,0]
      qt_fp = var(pm,nw, :qt_br)[f_idx,1]
      qt_tm = var(pm,nw, :qt_br)[t_idx,0]
      qt_tp = var(pm,nw, :qt_br)[t_idx,1]

      cref_qf=JuMP.@constraint(pm.model, q_fr - qd_fm + qd_fp + qt_fm - qt_fp == 0)
      cref_qt=JuMP.@constraint(pm.model, q_to - qd_tm + qd_tp + qt_tm - qt_tp == 0)
    end
    return cref_pf,cref_pt,cref_qf,cref_qt
end

function constraint_abs_branch_flow_ordering(pm::AbstractPowerModel, l::Int; nw::Int=pm.cnw)
    branch = ref(pm, nw, :branch, l)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (l, f_bus, t_bus)
    t_idx = (l, t_bus, f_bus)
    u_K = var(pm,nw,:u_K)
    u_ord_aux = var(pm,nw,:u_ord_aux,l)

    pt_fm = var(pm,nw, :pt_br)[f_idx,0]
    pt_tm = var(pm,nw, :pt_br)[t_idx,0]
    pt_fp = var(pm,nw, :pt_br)[f_idx,1]
    pt_tp = var(pm,nw, :pt_br)[t_idx,1]
    pd_fm = var(pm,nw, :pd_br)[f_idx,0]
    pd_tm = var(pm,nw, :pd_br)[t_idx,0]
    pd_fp = var(pm,nw, :pd_br)[f_idx,1]
    pd_tp = var(pm,nw, :pd_br)[t_idx,1]

    cref=nothing

    if haskey( var(pm,nw), :p ) && haskey( var(pm,nw), :q )
        qt_fm = var(pm,nw, :qt_br)[f_idx,0]
        qt_tm = var(pm,nw, :qt_br)[t_idx,0]
        qt_fp = var(pm,nw, :qt_br)[f_idx,1]
        qt_tp = var(pm,nw, :qt_br)[t_idx,1]
        qd_fm = var(pm,nw, :qd_br)[f_idx,0]
        qd_tm = var(pm,nw, :qd_br)[t_idx,0]
        qd_fp = var(pm,nw, :qd_br)[f_idx,1]
        qd_tp = var(pm,nw, :qd_br)[t_idx,1]
        cref=JuMP.@constraint(pm.model, -(pt_fm + pt_fp + pt_tm + pt_tp + qt_fm + qt_fp + qt_tm + qt_tp)
		    + (pd_fm + pd_fp + pd_tm + pd_tp + qd_fm + qd_fp + qd_tm + qd_tp) + u_ord_aux + u_K >= 0)
    elseif haskey( var(pm,nw), :p )
        cref=JuMP.@constraint(pm.model, -(pt_fm + pt_fp + pt_tm + pt_tp) + (pd_fm + pd_fp + pd_tm + pd_tp ) + u_ord_aux + u_K >= 0)
    end
    return cref
end


text:
module IncomeTaxCalculations

import Dates
import Dates: Date, now, TimeType, Year

import Model_Household: Person
import STBParameters: IncomeTaxSys
import TBComponents: TaxResult, calctaxdue, RateBands, delete_thresholds_up_to, *
import Parameters: @with_kw

using Definitions

export calc_income_tax, old_enough_for_mca, apply_allowance, ITResult
export calculate_company_car_charge

@with_kw mutable struct ITResult
    total_tax :: Real = 0.0
    taxable_income :: Real = 0.0
    adjusted_net_income :: Real = 0.0
    total_income :: Real = 0.0
    non_savings :: Real = 0.0
    allowance   :: Real = 0.0
    non_savings_band :: Integer = 0
    savings :: Real = 0.0
    savings_band :: Integer = 0
    dividends :: Real = 0.0
    dividend_band :: Integer = 0
    unused_allowance :: Real = 0.0
    mca :: Real = 0.0
    transferred_allowance :: Real = 0.0
end

## FIXME all these constants should ultimately be parameters
const MCA_DATE = Date(1935,4,6) # fixme make this a parameter

const SAVINGS_INCOME = Incomes_Dict(
    bank_interest => 1.0,
    bonds_and_gilts => 1.0,
    other_investment_income => 1.0
)

const DIVIDEND_INCOME = Incomes_Dict(
    stocks_shares => 1.0
)
const Exempt_Income = Incomes_Dict(
    individual_savings_account=>1.0,
    local_taxes=>1.0,
    free_school_meals => 1.0,
    dlaself_care => 1.0,
    dlamobility => 1.0,
    child_benefit => 1.0,
    pension_credit => 1.0,
    bereavement_allowance_or_widowed_parents_allowance_or_bereavement=> 1.0,
    armed_forces_compensation_scheme => 1.0, # FIXME not in my list check this
    war_widows_or_widowers_pension => 1.0,
    severe_disability_allowance => 1.0,
    attendence_allowance => 1.0,
    industrial_injury_disablement_benefit => 1.0,
    employment_and_support_allowance => 1.0,
    incapacity_benefit => 1.0,## taxable after 29 weeks,
    income_support => 1.0,
    maternity_allowance => 1.0,
    maternity_grant_from_social_fund => 1.0,
    funeral_grant_from_social_fund => 1.0,
    guardians_allowance => 1.0,
    winter_fuel_payments => 1.0,
    dwp_third_party_payments_is_or_pc => 1.0,
    dwp_third_party_payments_jsa_or_esa => 1.0,
    extended_hb => 1.0,
    working_tax_credit => 1.0,
    child_tax_credit => 1.0,
    working_tax_credit_lump_sum => 1.0,
    child_tax_credit_lump_sum => 1.0,
    housing_benefit => 1.0,
    universal_credit => 1.0,
    personal_independence_payment_daily_living => 1.0,
    personal_independence_payment_mobility => 1.0 )

function make_all_taxable()::Incomes_Dict
    eis = union(Set( keys( Exempt_Income )), Definitions.Expenses )
    all_t = Incomes_Dict()
    for i in instances(Incomes_Type)
        if ! (i ∈ eis )
            all_t[i]=1.0
        end
    end
    all_t
end

function calculate_company_car_charge(
    pers :: Person,
    calculator :: Function = guess_car_percentage_2020_21 ) :: Real
    value = max(0.0, pers.company_car_value-
        pers.company_car_contribution )
    if pers.fuel_supplied > 0.0
        value += sys.fuel_imputation
    end
    pct = calculator( pers.company_car_fuel_type )
    value * pct
end

function make_non_savings()::Incomes_Dict
    excl = union(Set(keys(DIVIDEND_INCOME)), Set( keys(SAVINGS_INCOME)))
    nsi = make_all_taxable()
    for i in excl
        delete!( nsi, i )
    end
    nsi
end

const NON_SAVINGS_INCOME = make_non_savings()
const ALL_TAXABLE = make_all_taxable()

"""
Very rough approximation to MCA age - ignores all months since we don't have that in a typical dataset
TODO maybe overload this with age as a Date?
"""
function old_enough_for_mca(
    age            :: Integer,
    model_run_date :: TimeType = now() ) :: Bool
    (model_run_date - Year(age)) < MCA_DATE
end

function calculate_allowance( pers::Person, sys :: IncomeTaxSys ) :: Real
    allowance = sys.personal_allowance
    if pers.registered_blind
        allowance += sys.blind_persons_allowance
    end
    allowance
end

function apply_allowance( allowance::Real, income::Real )::Tuple
    r = max( 0.0, income - allowance )
    allowance = max(0.0, allowance-income)
    allowance,r
end

"""

Complete(??) income tax calculation, based on the Scottish/UK 2019 system.
Mostly taken from Melville (2019) chs 2-4.

FIXME this is too long and needs broken up.

problems:

1. we do this in strict non-savings, savings, dividends order; see 2(11) for examples where it's now advantageous to use a different order
2.

returns a single total tax liabilty, plus multiple intermediate numbers
in the `intermediate` dict

"""
function calc_income_tax(
    pers   :: Person,
    sys    :: IncomeTaxSys,
    intermediate :: Dict,
    spouse_transfer :: Real = 0.0 ) :: ITResult
    itres :: ITResult = ITResult()
    total_income = ALL_TAXABLE*pers.income;
    non_savings = NON_SAVINGS_INCOME*pers.income;
    savings = SAVINGS_INCOME*pers.income;
    dividends = DIVIDEND_INCOME*pers.income;
    allowance = calculate_allowance( pers, sys )
    # allowance reductions goes here

    non_dividends = non_savings + savings

    adjusted_net_income = total_income
    adjusted_net_income += calculate_company_car_charge(pers)
    # ...

    non_savings_tax = TaxResult(0.0, 0)
    savings_tax = TaxResult(0.0, 0)
    dividend_tax = TaxResult(0.0, 0)

    if adjusted_net_income > sys.personal_allowance_income_limit
        allowance =
            max(0.0,
                allowance -
                    sys.personal_allowance_withdrawal_rate*(
                        adjusted_net_income - sys.personal_allowance_income_limit ))
    end
    taxable_income = adjusted_net_income-allowance
    intermediate["allowance"]=allowance
    intermediate["total_income"]=total_income
    intermediate["adjusted_net_income"]=adjusted_net_income
    intermediate["taxable_income"]=taxable_income
    intermediate["savings"]=savings
    intermediate["non_savings"]=non_savings
    intermediate["dividends"]=dividends

    savings_thresholds = deepcopy( sys.savings_thresholds )
    savings_rates = deepcopy( sys.savings_rates )
    # FIXME model all this with parameters
    toprate = size( savings_thresholds )[1]
    if taxable_income > 0
        allowance,non_savings_taxable = apply_allowance( allowance, non_savings )
        non_savings_tax = calctaxdue(
            taxable=non_savings_taxable,
            rates=sys.non_savings_rates,
            thresholds=sys.non_savings_thresholds )

        # horrific savings calculation see Melville Ch2 "Savings Income" & examples 2-3
        # FIXME Move to separate function
        # delete the starting bands up to non_savings taxabke icome
        savings_rates, savings_thresholds = delete_thresholds_up_to(
            rates=savings_rates,
            thresholds=savings_thresholds,
            upto=non_savings_taxable );
        if sys.personal_savings_allowance > 0
            psa = sys.personal_savings_allowance
            # println( "taxable income $taxable_income sys.savings_thresholds[2] $(sys.savings_thresholds[2])")
            if taxable_income > sys.savings_thresholds[toprate]
                psa = 0.0
            elseif taxable_income > sys.savings_thresholds[2] # above the basic rate
                psa *= 0.5 # FIXME parameterise this
            end
            if psa > 0.0 ## if we haven't deleted the zero band already, just widen it
                if savings_rates[1] == 0.0
                    savings_thresholds[1] += psa;
                else ## otherwise, insert a  new one.
                    savings_thresholds = vcat([psa], savings_thresholds )
                    savings_rates = vcat([0.0], savings_rates )
                end
            end
            intermediate["personal_savings_allowance"] = psa
        end # we have a personal_savings_allowance
        intermediate["savings_rates"] = savings_rates
        intermediate["savings_thresholds"] = savings_thresholds
        allowance,savings_taxable = apply_allowance( allowance, savings )
        savings_tax = calctaxdue(
            taxable=savings_taxable,
            rates=savings_rates,
            thresholds=savings_thresholds )

        # Dividends
        # see around example 8-9 ch2
        allowance,dividends_taxable =
            apply_allowance( allowance, dividends )
        dividend_rates=deepcopy(sys.dividend_rates)
        dividend_thresholds=deepcopy(sys.dividend_thresholds )
        # always preserve any bottom zero rate
        add_back_zero_band = false
        zero_band = 0.0
        used_thresholds = non_savings_taxable+savings_taxable
        copy_start = 1
        # handle the zero rate
        if dividend_rates[1] == 0.0
            add_back_zero_band = true
            zero_band = dividend_thresholds[1]
            used_thresholds += min( zero_band, dividends_taxable )
            copy_start = 2
        end
        dividend_rates, dividend_thresholds =
            delete_thresholds_up_to(
                rates=dividend_rates[copy_start:end],
                thresholds=dividend_thresholds[copy_start:end],
                upto=used_thresholds );
        if add_back_zero_band
            dividend_rates = vcat( [0.0], dividend_rates )
            dividend_thresholds .+= zero_band # push all up
            dividend_thresholds = vcat( zero_band, dividend_thresholds )
        end
        intermediate["dividend_rates"]=dividend_rates
        intermediate["dividend_thresholds"]=dividend_thresholds
        intermediate["add_back_zero_band"]=add_back_zero_band
        intermediate["dividends_taxable"]=dividends_taxable

        dividend_tax = calctaxdue(
            taxable=dividends_taxable,
            rates=dividend_rates,
            thresholds=dividend_thresholds )
    else # some allowance left
        allowance = -taxable_income # e.g. allowance - taxable_income
    end
    intermediate["non_savings_tax"]=non_savings_tax.due
    intermediate["savings_tax"]=savings_tax.due
    intermediate["dividend_tax"]=dividend_tax.due

    #
    # tax reducers
    #
    total_tax = non_savings_tax.due+savings_tax.due+dividend_tax.due
    if spouse_transfer > 0
        sp_reduction =
            sys.non_savings_rates[sys.non_savings_basic_rate]*spouse_transfer
        total_tax = max( 0.0, total_tax - sp_reduction )
    end
    itres.total_tax = total_tax
    itres.taxable_income = taxable_income
    itres.allowance = allowance
    itres.total_income = total_income
    itres.adjusted_net_income = adjusted_net_income
    itres.non_savings = non_savings_tax.due
    itres.non_savings_band = non_savings_tax.end_band
    itres.savings = savings_tax.due
    itres.savings_band = savings_tax.end_band
    itres.dividends = dividend_tax.due
    itres.dividend_band = dividend_tax.end_band
    itres.unused_allowance = allowance
    itres
end

function allowed_to_transfer_allowance(
    sys  :: IncomeTaxSys;
    from :: ITResult,
    to   :: ITResult ) :: Bool

   can_transfer :: Bool = true
   if ! (from.unused_allowance > 0.0 &&
         to.unused_allowance <= 0.0)
         # nothing to transfer - this is actually wrong since
         # you can opt to transfer some allowance even if you
         # can technically use it.
       can_transfer = false
   elseif to.savings_band > sys.savings_basic_rate ||
      to.non_savings_band > sys.non_savings_basic_rate ||
      to.dividend_band > sys.dividend_basic_rate
      can_transfer = false
   end
   ## TODO disallow if mca claimed
   can_transfer
end # can_transfer

function calculate_mca( pers :: Person, tax :: ITResult, sys :: IncomeTaxSys)::Real
    ## FIXME parameterise this
    mca = sys.married_couples_allowance
    if tax.adjusted_net_income > sys.mca_income_maximum
        mca = max( sys.mca_minimum, mca -
           (tax.adjusted_net_income-sys.mca_income_maximum)*sys.mca_withdrawal_rate)
    end
    mca * sys.mca_credit_rate
end

function calc_income_tax(
    head   :: Person,
    spouse :: Union{Nothing,Person},
    sys    :: IncomeTaxSys,
    intermediate :: Dict ) :: NamedTuple
    head_intermed = Dict()
    headtax = calc_income_tax( head, sys, head_intermed )
    intermediate["head_tax"] = head_intermed
    spousetax = nothing
    # FIXME the transferable stuff here
    # is not right as you can elect to transfer more than
    # the surplus allowance in some cases.
    # also - add in restrictions on transferring to
    # higher rate payers.
    if spouse != nothing
        spouse_intermed = Dict()
        intermediate["spouse_tax"] = spouse_intermed
        spousetax = calc_income_tax( spouse, sys, spouse_intermed )
        # This is not quite right - you can't claim the
        # MCA AND transfer an allowance. We're assuming
        # always MCA first (I think it's always more valuable?)
        if old_enough_for_mca( head.age ) || old_enough_for_mca( spouse.age )
            # shoud usually just go to the head but.. some stuff about partner
            # with greater income if married after 2005 and you can elect to do this if
            # married before, so:
            if headtax.adjusted_net_income > spousetax.adjusted_net_income
                headtax.mca = calculate_mca( head, headtax, sys )
                headtax.total_tax = max( 0.0, headtax.total_tax - headtax.mca )
            else
                spousetax.mca = calculate_mca( spouse, spousetax, sys )
                spousetax.total_tax = max( 0.0, spousetax.total_tax - spousetax.mca )
            end
        end
        if spousetax.mca == 0.0 == headtax.mca
            if allowed_to_transfer_allowance( sys, from=spousetax, to=headtax )
                transferable_allow = min( spousetax.unused_allowance, sys.marriage_allowance )
                headtax = calc_income_tax( head, sys, head_intermed, transferable_allow )
                intermediate["transfer_spouse_to_head"] = transferable_allow
            elseif allowed_to_transfer_allowance( sys, from=headtax, to=spousetax )
                transferable_allow = min( headtax.unused_allowance, sys.marriage_allowance )
                spousetax = calc_income_tax( spouse, sys, spouse_intermed, transferable_allow )
                intermediate["transfer_head_to_spouse"] = transferable_allow
            end
        end
    end
    ( head=headtax, spouse=spousetax )
end # calc_income_tax

end # module


text:
# This test is intended to be run after `setup.sh` runs. ExponentialUtilities
# should be held at v1.9.0.

using Revise, Pkg, Test

const thisdir = dirname(@__FILE__)
Pkg.activate(thisdir)
# This is only needed on Pkg versions that don't notify
Revise.active_project_watcher()

using ExponentialUtilities
id = Base.PkgId(ExponentialUtilities)
pkgdata = Revise.pkgdatas[id]
A = rand(3, 3); A = A'*A; A = A' + A;
@test_throws UndefVarError exponential!(A)   # not present on v1.9
# From a different process, switch the active version of ExponentialUtilities
run(Cmd(`julia -e 'using Pkg; Pkg.activate("."); Pkg.add(name="ExponentialUtilities", version="1.10.0")'`; dir=thisdir))
sleep(0.2)
revise()
@test exponential!(A) isa Matrix   # present on v1.9
# ...and then switch back (check that it's bidirectional and also to reset state)
run(Cmd(`julia -e 'using Pkg; Pkg.activate("."); Pkg.add(name="ExponentialUtilities", version="1.9.0")'`; dir=thisdir))
sleep(0.2)
revise()
@test_throws MethodError exponential!(A)   # not present on v1.9


text:
using MLJ, Tables, DataFrames, MLJDecisionTreeInterface, MLJLinearModels

models()

# 151-element Array{NamedTuple{(:name, :package_name, :is_supervised, :docstring, :hyperparameter_ranges, :hyperparameter_types, :hyperparameters, :implemented_methods, :is_pure_julia, :is_wrapper, :load_path, :package_license, :package_url, :package_uuid, :prediction_type, :supports_online, :supports_weights, :input_scitype, :target_scitype, :output_scitype),T} where T<:Tuple,1}:
#  (name = ARDRegressor, package_name = ScikitLearn, ... )
#  (name = AdaBoostClassifier, package_name = ScikitLearn, ... )
#  (name = AdaBoostRegressor, package_name = ScikitLearn, ... )
#  ⋮
#  (name = XGBoostClassifier, package_name = XGBoost, ... )
#  (name = XGBoostCount, package_name = XGBoost, ... )
#  (name = XGBoostRegressor, package_name = XGBoost, ... )

# use house price data from US Census Bureau
df = OpenML.load(574) |> DataFrame
X = df[:,[:P1,:P5p1,:P6p2,:P11p4,:P14p9,:P15p1,:P15p3,:P16p2,:P18p2,:P27p4,:H2p2,:H8p2,:H10p1,:H13p1,:H18pA,:H40p4]]
X = Tables.rowtable(X)
y = log.(df.price)

models(matching(X,y))

# declare a tree and lasso model
tree_model = @load DecisionTreeRegressor pkg=DecisionTree
lasso_model = @load LassoRegressor pkg=MLJLinearModels

# initialize "machines" where results can be reported
tree = machine(tree_model, X, y)
lass = machine(lasso_model, X, y)

# split into training and testing data
train, test = partition(eachindex(y), 0.7, shuffle=true)

# train the models
MLJ.fit!(tree, rows=train)
MLJ.fit!(lass, rows=train)

# predict in test set
yhat = MLJ.predict(tree, X[test,:]);
yhat = MLJ.predict(lass, X[test,:]);

# get RMSE across validation folds
MLJ.evaluate(tree_model,  X, y, resampling=CV(nfolds=6, shuffle=true), measure=rmse)
MLJ.evaluate(lasso_model, X, y, resampling=CV(nfolds=6, shuffle=true), measure=rmse)


text:
using Random
using Functors
using Zygote

include("grid.jl")
include("radfuncs.jl")
include("diffrules.jl")
Random.seed!(1)

mutable struct Op
    l
    kernel::AbstractArray
    grid::Grid
    radfunc
    rmin
    rmax
    boundary
end
@functor Op
Flux.trainable(m::Op) = [m.radfunc]

function makekernel(radfunc, rmin, rmax, l, grid)
    @unpack cell, coords, dv = grid
    n = size(cell, 1)
rnorms=0
Zygote.ignore() do
    rnorms = fieldnorm(coords)
end
    f = r-> rmin <= r <= rmax ? radfunc(r) : 0.0
    rscalars = f.(rnorms)
    if l == 0
        kernel = cat(rscalars, dims = n + 1)
    elseif l == 1
        kernel = broadcast(
            (a, r, s) -> (r == 0 || r > rmax) ? 0.0 : a * s / r,
            coords,
            repeat(rnorms, ones(Int, n)..., size(coords)[end]),
            repeat(rscalars, ones(Int, n)..., size(coords)[end]),
        )
    end
    dv * kernel
end
function Op(
    radfunc,
    rmin,
    rmax,
    cell::Matrix;
    l = 0,
    boundary =:zero,
)
    grid = Grid(cell, rmax)
    kernel = makekernel(radfunc, rmin, rmax, l, grid)

    Op(l, kernel, grid, radfunc,rmin, rmax, boundary)
end

"""
    Op(
        name::Union{Symbol,String},
        cell;
        boundary =:zero,
        rmin = 0,
        rmax = Inf,
        l = 0,
        σ = 1.0,
    )
    Op(
        radfunc,
        rmin::AbstractFloat,
        rmax::AbstractFloat,
        cell;
        l = 0,
        boundary =:zero,
    )

`Op` constructs finite difference operators. Prebuilt operators like differential operators (`▽`) & common Green's functions can be specified by name. Custom equivariant operators can be made by specifying radial function.
"""
function Op(
    name::Union{Symbol,String},
    cell;
    rmin = 0.,
    rmax = Inf,
    boundary =:zero,
    l = 0,
    σ = 1.0,
)
    name = Symbol(name)
    if name == :▽
        name = :grad
    end

    n = size(cell, 1)
    radfunc=nothing
    if name == :neural
    elseif name == :Gaussian
        radfunc = r -> exp(-r^2 / (2 * σ^2)) / sqrt(2π * σ^(2n))
        return Op(radfunc, rmin,3σ, cell)
    elseif name == :grad
        boundary = :smooth
        l = 1
        grid = Grid(cell, fill(3, n))

        if n == 1
            kernel = [1 / (2 * cell[1]), 0, -1 / (2 * cell[1])]
        elseif n == 2
            kernel = [
                sum(abs.([x, y])) > 1 ? zeros(2) : -cell' \ [x, y] / 2
                for x = -1:1, y = -1:1
            ]
            kernel = cat([getindex.(kernel, i) for i = 1:n]..., dims = n + 1)
        elseif n == 3
            kernel = [
                sum(abs.([x, y, z])) > 1 ? zeros(3) :
                    -cell' \ [x, y, z] / 2 for x = -1:1, y = -1:1, z = -1:1
            ]
            kernel = cat([getindex.(kernel, i) for i = 1:n]..., dims = n + 1)
        end
    end
    Op(l, kernel, grid, radfunc,rmin, rmax, boundary)
end

"""
    function (m::Op)(x::AbstractArray, )

"""
function (m::Op)(x::AbstractArray)
    li = 0
    lo = m.l
    m(x, li, lo)
end

function (m::Op)(x::AbstractArray, li, lo)
    @unpack grid, kernel, l, boundary = m
    @unpack origin, cell = grid
    ix = 0
    Zygote.ignore() do
        if boundary == :smooth
            ix = [
                Int.([a + 1, ((a+1):(b+a-2))..., b + a - 2])
                for (a, b) in zip(origin, size(x))
            ]
        else
            ix = [Int(a):Int(b + a - 1) for (a, b) in zip(origin, size(x))]
        end
    end

    FieldConv(li, l, lo)(x, kernel)[ix..., :]
end

function LinearAlgebra.:⋅(m::Op, x)
    m(x, m.l, 0)
end
function ⨉(m::Op, x)
    m(x, 1, 1)
end
function remake!(m)
    @unpack radfunc,rmin, rmax, l, grid = m
    m.kernel = makekernel(radfunc,rmin, rmax, l, grid)
end


text:
# returns all variables of a cost function, in terms of appearance
extract_variables(t::Term) = variables(t)

function extract_variables(t::NTuple{N,Term}) where {N}
	x = variables.(t)
	xAll = x[1]
	for i = 2:length(x)
		for xi in x[i]
			if (xi in xAll) == false
				xAll = (xAll...,xi)
			end
		end
	end
	return xAll
end

# extract functions from terms
function extract_functions(t::Term)
	f = displacement(t) == 0 ? t.f : PrecomposeDiagonal(t.f, 1.0, displacement(t)) #for now I keep this
	f = t.lambda == 1. ? f : Postcompose(f, t.lambda)                                  #for now I keep this
	#TODO change this
	return f
end
extract_functions(t::NTuple{N,Term}) where {N} = SeparableSum(extract_functions.(t))
extract_functions(t::Tuple{Term}) = extract_functions(t[1])

# extract functions from terms without displacement
function extract_functions_nodisp(t::Term)
	f = t.lambda == 1. ? t.f : Postcompose(t.f, t.lambda)
	return f
end
extract_functions_nodisp(t::NTuple{N,Term}) where {N} = SeparableSum(extract_functions_nodisp.(t))
extract_functions_nodisp(t::Tuple{Term}) = extract_functions_nodisp(t[1])

# extract operators from terms

# returns all operators with an order dictated by xAll

#single term, single variable
extract_operators(xAll::Tuple{Variable}, t::Term)  = operator(t)

extract_operators(xAll::NTuple{N,Variable}, t::Term) where {N} = extract_operators(xAll, (t,))

#multiple terms, multiple variables
function extract_operators(xAll::NTuple{N,Variable}, t::NTuple{M,Term}) where {N,M}
	ops = ()
	for ti in t
		tex = expand(xAll,ti)
		ops = (ops...,sort_and_extract_operators(xAll,tex))
	end
	return vcat(ops...)
end

sort_and_extract_operators(xAll::Tuple{Variable}, t::Term) = operator(t)

function sort_and_extract_operators(xAll::NTuple{N,Variable}, t::Term) where {N}
	p = zeros(Int,N)
	xL = variables(t)
	for i in eachindex(xAll)
		p[i] = findfirst( xi -> xi == xAll[i], xL)
	end
	return operator(t)[p]
end

# extract affines from terms

# returns all affines with an order dictated by xAll

#single term, single variable
extract_affines(xAll::Tuple{Variable}, t::Term)  = affine(t)

extract_affines(xAll::NTuple{N,Variable}, t::Term) where {N} = extract_affines(xAll, (t,))

#multiple terms, multiple variables
function extract_affines(xAll::NTuple{N,Variable}, t::NTuple{M,Term}) where {N,M}
	ops = ()
	for ti in t
		tex = expand(xAll,ti)
		ops = (ops...,sort_and_extract_affines(xAll,tex))
	end
	return vcat(ops...)
end

sort_and_extract_affines(xAll::Tuple{Variable}, t::Term) = affine(t)

function sort_and_extract_affines(xAll::NTuple{N,Variable}, t::Term) where {N}
	p = zeros(Int,N)
	xL = variables(t)
	for i in eachindex(xAll)
		p[i] = findfirst( xi -> xi == xAll[i], xL)
	end
	return affine(t)[p]
end

# expand term domain dimensions
function expand(xAll::NTuple{N,Variable}, t::Term{T1,T2,T3}) where {N,T1,T2,T3}
	xt   = variables(t)
	C    = codomainType(operator(t))
	size_out = size(operator(t),1)
	ex = t.A

	for x in xAll
		if !( x in variables(t) )
			ex += Zeros(eltype(~x),size(x),C,size_out)*x
		end
	end
	return Term(t.lambda, t.f, ex)
end

# extract function and merge operator
function extract_merge_functions(t::Term)
    if is_sliced(t)
        if typeof(operator(t)) <: Compose
            op = operator(t).A[2]
        else
            op = Eye(size(operator(t),1)...)
        end
    else
        op = operator(t)
    end
    if is_eye(op)
        f = displacement(t) == 0 ? t.f : PrecomposeDiagonal(t.f, 1.0, displacement(t))
    elseif is_diagonal(op)
        f = PrecomposeDiagonal(t.f, diag(op), displacement(t))
    elseif is_AAc_diagonal(op)
        f = Precompose(t.f, op, diag_AAc(op), displacement(t))
    end
	f = t.lambda == 1. ? f : Postcompose(f, t.lambda) #for now I keep this
	#TODO change this
	return f
end

function extract_proximable(xAll::NTuple{N,Variable}, t::NTuple{M,Term}) where {N,M}
	fs = ()
	for x in xAll
		tx = () #terms containing x
		for ti in t
			if x in variables(ti)
				tx = (tx...,ti) #collect terms containing x
			end
		end
		if isempty(tx)
			fx = IndFree()
		elseif length(tx) == 1          #only one term per variable
			fx = extract_proximable(x,tx[1])
		else
            #multiple terms per variable
            #currently this happens only with GetIndex
			fxi,idxs = (),()
			for ti in tx
				fxi  = (fxi..., extract_merge_functions(ti))
                idx = typeof(operator(ti)) <: Compose ? operator(ti).A[1].idx : operator(ti).idx
				idxs = (idxs...,  idx   )
			end
			fx = SlicedSeparableSum(fxi,idxs)
		end
		fs = (fs...,fx)
	end
	if length(fs) > 1
		return SeparableSum(fs)  ##probably change constructor in Prox?
	else
		return fs[1]
	end
end

extract_proximable(xAll::Variable, t::Term) =  extract_merge_functions(t)
extract_proximable(xAll::NTuple{N,Variable}, t::Term) where {N} = extract_proximable(xAll,(t,))


text:
## Author:  Pearson
## Date: April 2020
## Julia version: 1.3.1
## Purpose: Solve EKNR(2016) Full Model - Baseline Equilibrium

# Precompile Packages
using Parameters
using Tables
using CSV
using NLsolve
#using Plots
include("Full_Baseline_Functions.jl")
## Read data
Data = CSV.read(joinpath(@__DIR__, "..", "input", "part1_21cty.csv"))
π_data = CSV.read(joinpath(@__DIR__, "..", "input", "part2_21cty.csv"))

###################################################################
# Parameters
###################################################################
# Program Parameters
NC = 21 # Number of countries
NS = 3 # Number of sectors
NK = 2 # Number of capital types(C and D, corresponding to sector 1 and 2)
T₁ = 211
T_data = 0
T_tail = 150
T = T_data+T_tail

# Pick the final period
Data = Data[Data.date.==T₁, 3:end]
π_data = π_data[π_data.date.==T₁, 5:end]

# Preference
ρ = Data.rho[1] # Intertemporal elasticity of substitution
ψ = [Data.psiC[1] Data.psiD[1] Data.psiN[1]]

# Production
α = [Data.alphaC[1] Data.alphaD[1] 0.5] # Investment efficiency
δ = [Data.deltaC[1] Data.deltaD[1]] # Depreciation rete
θ = 2 # Trade elasticity for sector D and S
β̃ᴸ = ones(NC, NS+1)
β̃ᴷ = ones(NC, NS+1, NK)
β̃ᴹ = ones(NC, NS+1, NS)
β̃ᴸ = [Data.betaTilLC Data.betaTilLD Data.betaTilLS Data.betaTilLR] # Labor's income share
β̃ᴷ_temp = [Data.betaTilKCC Data.betaTilKDC Data.betaTilKSC Data.betaTilKRC Data.betaTilKCD Data.betaTilKDD Data.betaTilKSD Data.betaTilKRD]
β̃ᴹ_temp = [Data.betaTilICC Data.betaTilIDC Data.betaTilISC Data.betaTilIRC Data.betaTilICD Data.betaTilIDD Data.betaTilISD Data.betaTilIRD Data.betaTilICS Data.betaTilIDS Data.betaTilISS Data.betaTilIRS]
for country in eachindex(β̃ᴷ_temp[:,1])
    β̃ᴷ[country, :, :] = reshape(β̃ᴷ_temp[country, :], NS+1, NK) # Capital's income share
    β̃ᴹ[country, :, :] = reshape(β̃ᴹ_temp[country, :], NS+1, NS) # Intermediate's income share
end
###################################################################
# Guess
###################################################################
K̂ = ones(NC, NK, T)
Ŷ = ones(NC, NK, T)
guess_dynamic = Array{Float64,3}(undef,NC,NK,T)
guess_dynamic[1:NC,1:NK,1] = K̂[1:NC,1:NK,1]
guess_dynamic[1:NC,1:NK,2:T] = Ŷ[1:NC,1:NK,1:T-1]
guess_1d = reshape(guess_dynamic,NC*NK*T)
###################################################################
# Preallocate memory
###################################################################
# Level variables
π = zeros(NC, NC, NS, T)
Y = zeros(NC, NS+1, T)
# X = zeros(NC, NS+1, T)
Xᶠ = zeros(NC, NS+1, T)
rK = zeros(NC, NK, T)
wL = zeros(NC, T)

# Exogenous variables
D = zeros(NC, NS+1)

###################################################################
# Initial Conditions
###################################################################
for country in 1:NC
    π[1:NC,country,1:NS,1] = Matrix(π_data[(country-1)*NC+1:country*NC, 2:4])
end
Y[1:NC,1:NS+1,1] = [Data.yC1 Data.yD1 Data.yS1 Data.yR1]
# X[1:NC,1:NS+1,1] = [Data.xC1 Data.xD1 Data.xS1 Data.xR1]
Xᶠ[1:NC,1:NS+1,1] = [Data.xFC1 Data.xFD1 Data.xFS1 Data.xFR1]
rK[1:NC,1:NK,1] = [Data.capCinc1 Data.capDinc1]
wL[1:NC,1] = Data.labinc1
X = zeros(NC,NS,1)
for s in 1:NS
    X[1:NC,s,1] = π[1:NC,1:NC,s,1]'\(Y[1:NC,s,1])
end
###################################################################
# Exogenous Variables
###################################################################
# Time-varying variables(Hold constant in this procedure)
Â = ones(NC, NS, T)
χ̂ = ones(NC, NS, T)
d̂ = ones(NC, NC, NS, T)
L̂ = ones(NC, T)
ϕ̂ = ones(NC, T)

D[1:NC,1:NS+1,1] = [Data.DC1 Data.DD1 Data.DS1 Data.DR1]
D = repeat(D,1,1,T)
T̂ = ones(NC,NS,T)
#ψ̂ˢ = ones(NC, T)
#ϕ̂ᵂ = ones(NC*NS, T)

## Pack init conditions, exogenous variables and parameters
myinit_dynamic = @with_kw (
    π₁ = π[1:NC,1:NC,1:NS,1], Y₁ = Y[1:NC,1:NS+1,1], Xᶠ₁ = Xᶠ[1:NC,1:NS+1,1], wL₁ = wL[1:NC,1], rK₁ = rK[1:NC,1:NK,1])
myexos_dynamic = @with_kw (
    Dᴿ = D[1:NC,4,1:T], L̂ = L̂[1:NC,1:T], d̂ = d̂[1:NC,1:NC,1:NS,1:T], T̂ = T̂[1:NC,1:NS,1:T])
myparams_dynamic = @with_kw (
    T = T, NC = NC, NS = NS, NK = NK, β̃ᴸ = β̃ᴸ, β̃ᴷ = β̃ᴷ, ψ = ψ, θ = θ, β̃ᴹ = β̃ᴹ, ρ = ρ, δ = δ, α = α)
init_dynamic = myinit_dynamic()
exos_dynamic = myexos_dynamic()
params_dynamic = myparams_dynamic()
## Solve the dynamic problem
# Solve the fix point problem
println("Start to solve the dynamic problem.")
println("Run time and memory cost:")
include("Full_Baseline_Functions.jl")
res_dynamic = similar(guess_dynamic)
@time results_dynamic =
    try
        results_dynamic = nlsolve(
            (res_dynamic, guess_dynamic) -> dynamic_problem!(
                res_dynamic, guess_dynamic, init_dynamic, exos_dynamic, params_dynamic),
            guess_1d,
            ftol=1e-6,
            method=:newton,
            autodiff=:forward,
            show_trace=true,
        )
    catch err
        if isa(err, DomainError)
            error("Failed to solve the dynamic problem, please try again.")
        end
    end

# Check Convergence
converged(results_dynamic) || error("Failed to converge in $(results_dynamic.iterations) iterations.")
println("Successfully solved the dynamic problem.\n")

# Catch Solutions
res_dynamic = similar(results_dynamic.zero)
res_dynamic, K̂, Ŷ = dynamic_problem!(
    res_dynamic, results_dynamic.zero, exos_dynamic, params_dynamic)

## Config use
show(Data[1,:], allcols=:true)


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

using Test

import Base: root_module

Foo_module = :Foo4b3a94a1a081a8cb
Foo2_module = :F2oo4b3a94a1a081a8cb
FooBase_module = :FooBase4b3a94a1a081a8cb
@eval module ConflictingBindings
    export $Foo_module, $FooBase_module
    $Foo_module = 232
    $FooBase_module = 9134
end
using .ConflictingBindings

# this environment variable would affect some error messages being tested below
# so we disable it for the tests below
withenv( "JULIA_DEBUG_LOADING" => nothing ) do

dir = mktempdir()
dir2 = mktempdir()
insert!(LOAD_PATH, 1, dir)
insert!(Base.LOAD_CACHE_PATH, 1, dir)
try
    Foo_file = joinpath(dir, "$Foo_module.jl")
    Foo2_file = joinpath(dir, "$Foo2_module.jl")
    FooBase_file = joinpath(dir, "$FooBase_module.jl")

    write(FooBase_file,
          """
          __precompile__(true)

          module $FooBase_module
              import Base: hash, >
              struct fmpz end
              struct typeA end
              >(x::fmpz, y::Int) = Base.cmp(x, y) > 0
              function hash(a::typeA, h::UInt)
                  d = den(a)
                  return h
              end
          end
          """)
    write(Foo2_file,
          """
          __precompile__(true)

          module $Foo2_module
              export override
              override(x::Integer) = 2
              override(x::AbstractFloat) = Float64(override(1))
          end
          """)
    write(Foo_file,
          """
          __precompile__(true)

          module $Foo_module
              using $FooBase_module, $FooBase_module.typeA
              import $Foo2_module: $Foo2_module, override
              import $FooBase_module.hash
              import Test

              struct typeB
                  y::typeA
              end
              hash(x::typeB) = hash(x.y)

              # test that docs get reconnected
              @doc "foo function" foo(x) = x + 1
              include_dependency("foo.jl")
              include_dependency("foo.jl")
              module Bar
                  @doc "bar function" bar(x) = x + 2
                  include_dependency("bar.jl")
              end

              # test for creation of some reasonably complicated type
              struct MyType{T} end
              const t17809s = Any[
                    Tuple{
                        Type{Ptr{MyType{i}}},
                        Array{Ptr{MyType{MyType{:sym}()}}(0), 0},
                        Val{Complex{Int}(1, 2)},
                        Val{3},
                        Val{nothing}}
                    for i = 0:25]

              # test that types and methods get reconnected correctly
              # issue 16529 (adding a method to a type with no instances)
              (::Task)(::UInt8, ::UInt16, ::UInt32) = 2

              # issue 16471 (capturing references to a kwfunc)
              Test.@test_throws ErrorException Core.kwfunc(Base.nothing)
              Base.nothing(::UInt8, ::UInt16, ::UInt32; x = 52) = x
              const nothingkw = Core.kwfunc(Base.nothing)

              # issue 16908 (some complicated types and external method definitions)
              abstract type CategoricalPool{T, R <: Integer, V} end
              abstract type CategoricalValue{T, R <: Integer} end
              struct NominalPool{T, R <: Integer, V} <: CategoricalPool{T, R, V}
                  index::Vector{T}
                  invindex::Dict{T, R}
                  order::Vector{R}
                  ordered::Vector{T}
                  valindex::Vector{V}
              end
              struct NominalValue{T, R <: Integer} <: CategoricalValue{T, R}
                  level::R
                  pool::NominalPool{T, R, NominalValue{T, R}}
              end
              struct OrdinalValue{T, R <: Integer} <: CategoricalValue{T, R}
                  level::R
                  pool::NominalPool{T, R, NominalValue{T, R}}
              end
              (::Union{Type{NominalValue}, Type{OrdinalValue}})() = 1
              (::Union{Type{NominalValue{T}}, Type{OrdinalValue{T}}})() where {T} = 2
              (::Type{Vector{NominalValue{T, R}}})() where {T, R} = 3
              (::Type{Vector{NominalValue{T, T}}})() where {T} = 4
              (::Type{Vector{NominalValue{Int, Int}}})() = 5

              # more tests for method signature involving a complicated type
              # issue 18343
              struct Pool18343{R, V}
                  valindex::Vector{V}
              end
              struct Value18343{T, R}
                  pool::Pool18343{R, Value18343{T, R}}
              end
              Base.convert(::Type{Nullable{S}}, ::Value18343{Nullable}) where {S} = 2
              Base.convert(::Type{Nullable{Value18343}}, ::Value18343{Nullable}) = 2
              Base.convert(::Type{Ref}, ::Value18343{T}) where {T} = 3

              let some_method = @which Base.include("string")
                    # global const some_method // FIXME: support for serializing a direct reference to an external Method not implemented
                  global const some_linfo =
                      ccall(:jl_specializations_get_linfo, Ref{Core.MethodInstance}, (Any, Any, Any, UInt),
                          some_method, Tuple{typeof(Base.include), String}, Core.svec(), typemax(UInt))
              end

              g() = override(1.0)
              Test.@test g() === 2.0 # compile this
          end
          """)
    @test_throws ErrorException Core.kwfunc(Base.nothing) # make sure `nothing` didn't have a kwfunc (which would invalidate the attempted test)

    # Issue #12623
    @test __precompile__(true) === nothing

    # Issue #21307
    Base.require(Foo2_module)
    @eval let Foo2_module = $(QuoteNode(Foo2_module)), # use @eval to see the results of loading the compile
              Foo = root_module(Foo2_module)
        Foo.override(::Int) = 'a'
        Foo.override(::Float32) = 'b'
    end

    Base.require(Foo_module)

    @eval let Foo_module = $(QuoteNode(Foo_module)), # use @eval to see the results of loading the compile
              Foo = root_module(Foo_module)
        @test Foo.foo(17) == 18
        @test Foo.Bar.bar(17) == 19

        # Issue #21307
        @test Foo.g() === 97.0
        @test Foo.override(1.0e0) == Float64('a')
        @test Foo.override(1.0f0) == 'b'
        @test Foo.override(UInt(1)) == 2
    end

    cachefile = joinpath(dir, "$Foo_module.ji")
    # use _require_from_serialized to ensure that the test fails if
    # the module doesn't reload from the image:
    @test_warn "WARNING: replacing module $Foo_module." begin
        ms = Base._require_from_serialized(Foo_module, cachefile)
        @test isa(ms, Array{Any,1})
        Base.register_all(ms)
    end

    let Foo = root_module(Foo_module)
        @test_throws MethodError Foo.foo(17) # world shouldn't be visible yet
    end
    @eval let Foo_module = $(QuoteNode(Foo_module)), # use @eval to see the results of loading the compile
              Foo2_module = $(QuoteNode(Foo2_module)),
              FooBase_module = $(QuoteNode(FooBase_module)),
              Foo = root_module(Foo_module),
              dir = $(QuoteNode(dir)),
              cachefile = $(QuoteNode(cachefile)),
              Foo_file = $(QuoteNode(Foo_file))
        @test Foo.foo(17) == 18
        @test Foo.Bar.bar(17) == 19

        # Issue #21307
        @test Foo.g() === 97.0
        @test Foo.override(1.0e0) == Float64('a')
        @test Foo.override(1.0f0) == 'b'
        @test Foo.override(UInt(1)) == 2

        # issue #12284:
        @test stringmime("text/plain", Base.Docs.doc(Foo.foo)) == "foo function\n"
        @test stringmime("text/plain", Base.Docs.doc(Foo.Bar.bar)) == "bar function\n"

        modules, deps, required_modules = Base.parse_cache_header(cachefile)
        discard_module = mod_fl_mt -> (mod_fl_mt[2], mod_fl_mt[3])
        @test modules == Dict(Foo_module => Base.module_uuid(Foo))
        @test map(x -> x[1],  sort(discard_module.(deps))) == [Foo_file, joinpath(dir, "bar.jl"), joinpath(dir, "foo.jl")]
        srctxt = Base.read_dependency_src(cachefile, Foo_file)
        @test !isempty(srctxt) && srctxt == read(Foo_file, String)
        @test_throws ErrorException Base.read_dependency_src(cachefile, "/tmp/nonexistent.txt")
        # dependencies declared with `include_dependency` should not be stored
        @test_throws ErrorException Base.read_dependency_src(cachefile, joinpath(dir, "foo.jl"))

        modules, deps1 = Base.cache_dependencies(cachefile)
        @test modules == merge(Dict(s => Base.module_uuid(getfield(Foo, s)) for s in
                                    [:Base, :Core, Foo2_module, FooBase_module, :Main, :Test]),
                               # plus modules included in the system image
                               Dict(s => Base.module_uuid(Base.root_module(s)) for s in
                                    [:DelimitedFiles,:Mmap]))
        @test discard_module.(deps) == deps1

        @test current_task()(0x01, 0x4000, 0x30031234) == 2
        @test nothing(0x01, 0x4000, 0x30031234) == 52
        @test nothing(0x01, 0x4000, 0x30031234; x = 9142) == 9142
        @test Foo.nothingkw === Core.kwfunc(Base.nothing)

        @test Foo.NominalValue() == 1
        @test Foo.OrdinalValue() == 1
        @test Foo.NominalValue{Int}() == 2
        @test Foo.OrdinalValue{Int}() == 2
        let T = Vector{Foo.NominalValue{Int}}
            @test isa(T(), T)
        end
        @test Vector{Foo.NominalValue{Int32, Int64}}() == 3
        @test Vector{Foo.NominalValue{UInt, UInt}}() == 4
        @test Vector{Foo.NominalValue{Int, Int}}() == 5
        @test all(i -> Foo.t17809s[i + 1] ===
            Tuple{
                Type{Ptr{Foo.MyType{i}}},
                Array{Ptr{Foo.MyType{Foo.MyType{:sym}()}}(0), 0},
                Val{Complex{Int}(1, 2)},
                Val{3},
                Val{nothing}},
            0:25)
        some_method = @which Base.include("string")
        some_linfo =
                ccall(:jl_specializations_get_linfo, Ref{Core.MethodInstance}, (Any, Any, Any, UInt),
                    some_method, Tuple{typeof(Base.include), String}, Core.svec(), typemax(UInt))
        @test Foo.some_linfo::Core.MethodInstance === some_linfo

        PV = Foo.Value18343{Nullable}.body.types[1]
        VR = PV.types[1].parameters[1]
        @test PV.types[1] === Array{VR,1}
        @test pointer_from_objref(PV.types[1]) ===
              pointer_from_objref(PV.types[1].parameters[1].types[1].types[1])
        @test PV === PV.types[1].parameters[1].types[1]
        @test pointer_from_objref(PV) === pointer_from_objref(PV.types[1].parameters[1].types[1])
    end

    Baz_file = joinpath(dir, "Baz.jl")
    write(Baz_file,
          """
          __precompile__(false)
          module Baz
          end
          """)

    @test_warn "ERROR: LoadError: Declaring __precompile__(false) is not allowed in files that are being precompiled.\nStacktrace:\n [1] __precompile__" try
        Base.compilecache("Baz") # from __precompile__(false)
        error("__precompile__ disabled test failed")
    catch exc
        isa(exc, ErrorException) || rethrow(exc)
        !isempty(search(exc.msg, "__precompile__(false)")) && rethrow(exc)
    end

    # Issue #12720
    FooBar1_file = joinpath(dir, "FooBar1.jl")
    write(FooBar1_file,
          """
          __precompile__(true)
          module FooBar1
              using FooBar
          end
          """)
    sleep(2) # give FooBar and FooBar1 different timestamps, in reverse order too
    FooBar_file = joinpath(dir, "FooBar.jl")
    write(FooBar_file,
          """
          __precompile__(true)
          module FooBar
          end
          """)

    Base.compilecache("FooBar")
    @test isfile(joinpath(dir, "FooBar.ji"))
    @test !Base.stale_cachefile(FooBar_file, joinpath(dir, "FooBar.ji"))
    @test !isdefined(Main, :FooBar)
    @test !isdefined(Main, :FooBar1)

    relFooBar_file = joinpath(dir, "subfolder", "..", "FooBar.jl")
    @test Base.stale_cachefile(relFooBar_file, joinpath(dir, "FooBar.ji")) == !Sys.iswindows() # `..` is not a symlink on Windows
    mkdir(joinpath(dir, "subfolder"))
    @test !Base.stale_cachefile(relFooBar_file, joinpath(dir, "FooBar.ji"))

    @eval using FooBar
    fb_uuid = Base.module_uuid(FooBar)
    sleep(2); touch(FooBar_file)
    insert!(Base.LOAD_CACHE_PATH, 1, dir2)
    @test Base.stale_cachefile(FooBar_file, joinpath(dir, "FooBar.ji"))
    @eval using FooBar1
    @test !isfile(joinpath(dir2, "FooBar.ji"))
    @test !isfile(joinpath(dir, "FooBar1.ji"))
    @test isfile(joinpath(dir2, "FooBar1.ji"))
    @test Base.stale_cachefile(FooBar_file, joinpath(dir, "FooBar.ji"))
    @test !Base.stale_cachefile(FooBar1_file, joinpath(dir2, "FooBar1.ji"))
    @test fb_uuid == Base.module_uuid(FooBar)
    fb_uuid1 = Base.module_uuid(FooBar1)
    @test fb_uuid != fb_uuid1

    reload("FooBar")
    @test fb_uuid != Base.module_uuid(root_module(:FooBar))
    @test fb_uuid1 == Base.module_uuid(FooBar1)
    fb_uuid = Base.module_uuid(root_module(:FooBar))
    @test isfile(joinpath(dir2, "FooBar.ji"))
    @test Base.stale_cachefile(FooBar_file, joinpath(dir, "FooBar.ji"))
    @test !Base.stale_cachefile(FooBar1_file, joinpath(dir2, "FooBar1.ji"))
    @test !Base.stale_cachefile(FooBar_file, joinpath(dir2, "FooBar.ji"))

    reload("FooBar1")
    @test fb_uuid == Base.module_uuid(root_module(:FooBar))
    @test fb_uuid1 != Base.module_uuid(root_module(:FooBar1))

    @test isfile(joinpath(dir2, "FooBar.ji"))
    @test isfile(joinpath(dir2, "FooBar1.ji"))
    @test Base.stale_cachefile(FooBar_file, joinpath(dir, "FooBar.ji"))
    @test !Base.stale_cachefile(FooBar_file, joinpath(dir2, "FooBar.ji"))
    @test !Base.stale_cachefile(FooBar1_file, joinpath(dir2, "FooBar1.ji"))

    # test checksum
    open(joinpath(dir2, "FooBar1.ji"), "a") do f
        write(f, 0x076cac96) # append 4 random bytes
    end
    @test Base.stale_cachefile(FooBar1_file, joinpath(dir2, "FooBar1.ji"))

    # test behavior of precompile modules that throw errors
    FooBar2_file = joinpath(dir, "FooBar2.jl")
    write(FooBar2_file,
          """
          __precompile__(true)
          module FooBar2
          error("break me")
          end
          """)
    @test_warn "ERROR: LoadError: break me\nStacktrace:\n [1] error" try
        Base.require(:FooBar2)
        error("\"LoadError: break me\" test failed")
    catch exc
        isa(exc, ErrorException) || rethrow(exc)
        !isempty(search(exc.msg, "ERROR: LoadError: break me")) && rethrow(exc)
    end

    # Test transitive dependency for #21266
    FooBarT_file = joinpath(dir, "FooBarT.jl")
    write(FooBarT_file,
          """
          __precompile__(true)
          module FooBarT
          end
          """)
    FooBarT1_file = joinpath(dir, "FooBarT1.jl")
    write(FooBarT1_file,
          """
          __precompile__(true)
          module FooBarT1
              using FooBarT
          end
          """)
    FooBarT2_file = joinpath(dir, "FooBarT2.jl")
    write(FooBarT2_file,
          """
          __precompile__(true)
          module FooBarT2
              using FooBarT1
          end
          """)
    Base.compilecache("FooBarT2")
    write(FooBarT1_file,
          """
          __precompile__(true)
          module FooBarT1
          end
          """)
    rm(FooBarT_file)
    @test Base.stale_cachefile(FooBarT2_file, joinpath(dir2, "FooBarT2.ji"))
    @test Base.require(:FooBarT2) isa Module
finally
    splice!(Base.LOAD_CACHE_PATH, 1:2)
    splice!(LOAD_PATH, 1)
    rm(dir, recursive=true)
    rm(dir2, recursive=true)
end

# test --compiled-modules=no command line option
let dir = mktempdir(),
    Time_module = :Time4b3a94a1a081a8cb

    try
        write(joinpath(dir, "$Time_module.jl"),
              """
              module $Time_module
                  __precompile__(true)
                  time = Base.time()
              end
              """)

        eval(quote
            insert!(LOAD_PATH, 1, $(dir))
            insert!(Base.LOAD_CACHE_PATH, 1, $(dir))
            Base.compilecache(:Time4b3a94a1a081a8cb)
        end)

        exename = `$(Base.julia_cmd()) --compiled-modules=yes --startup-file=no`

        testcode = """
            insert!(LOAD_PATH, 1, $(repr(dir)))
            insert!(Base.LOAD_CACHE_PATH, 1, $(repr(dir)))
            using $Time_module
            getfield($Time_module, :time)
        """

        t1_yes = readchomp(`$exename --compiled-modules=yes -E $(testcode)`)
        t2_yes = readchomp(`$exename --compiled-modules=yes -E $(testcode)`)
        @test t1_yes == t2_yes

        t1_no = readchomp(`$exename --compiled-modules=no -E $(testcode)`)
        t2_no = readchomp(`$exename --compiled-modules=no -E $(testcode)`)
        @test t1_no != t2_no
        @test parse(Float64, t1_no) < parse(Float64, t2_no)

    finally
        splice!(Base.LOAD_CACHE_PATH, 1)
        splice!(LOAD_PATH, 1)
        rm(dir, recursive=true)
    end
end

# test loading a package with conflicting namespace
let dir = mktempdir()
    Test_module = :Test6c92f26
    try
        write(joinpath(dir, "Iterators.jl"),
              """
              module Iterators
                   __precompile__(true)
              end
              """)

        write(joinpath(dir, "$Test_module.jl"),
              """
              module $Test_module
                   __precompile__(true)
                   using Iterators
              end
              """)

        testcode = """
            insert!(LOAD_PATH, 1, $(repr(dir)))
            insert!(Base.LOAD_CACHE_PATH, 1, $(repr(dir)))
            using $Test_module
        """

        exename = `$(Base.julia_cmd()) --startup-file=no`
        let fname = tempname()
            try
                @test readchomp(pipeline(`$exename -E $(testcode)`, stderr=fname)) == "nothing"
                @test Test.ismatch_warn("WARNING: replacing module $Test_module.\n", read(fname, String))
            finally
                rm(fname, force=true)
            end
        end
        # Loading $Test_module from the cache should not bring `Base.Iterators`
        # into `Main`, since that would lead to a namespace conflict with
        # the module `Iterators` defined above.
        let fname = tempname()
            try
                @test readchomp(pipeline(`$exename -E $(testcode)`, stderr=fname)) == "nothing"
                # e.g `@test_nowarn`
                @test Test.ismatch_warn(r"^(?!.)"s, read(fname, String))
            finally
                rm(fname, force=true)
            end
        end
    finally
        rm(dir, recursive=true)
    end
end

let dir = mktempdir()
    try
        insert!(LOAD_PATH, 1, dir)
        insert!(Base.LOAD_CACHE_PATH, 1, dir)

        loaded_modules = Channel{Symbol}(32)
        callback = (mod::Symbol) -> put!(loaded_modules, mod)
        push!(Base.package_callbacks, callback)

        Test1_module = :Teste4095a81
        Test2_module = :Teste4095a82
        Test3_module = :Teste4095a83

        write(joinpath(dir, "$(Test1_module).jl"),
              """
              module $(Test1_module)
                  __precompile__(true)
              end
              """)

        Base.compilecache("$(Test1_module)")
        write(joinpath(dir, "$(Test2_module).jl"),
              """
              module $(Test2_module)
                  __precompile__(true)
                  using $(Test1_module)
              end
              """)
        Base.compilecache("$(Test2_module)")
        @test !Base.isbindingresolved(Main, Test2_module)
        Base.require(Test2_module)
        @test take!(loaded_modules) == Test1_module
        @test take!(loaded_modules) == Test2_module
        write(joinpath(dir, "$(Test3_module).jl"),
              """
              module $(Test3_module)
                  using $(Test3_module)
              end
              """)
        Base.require(Test3_module)
        @test take!(loaded_modules) == Test3_module
    finally
        pop!(Base.package_callbacks)
        splice!(Base.LOAD_CACHE_PATH, 1)
        splice!(LOAD_PATH, 1)
        rm(dir, recursive=true)
    end
end

let module_name = string("a",randstring())
    mktempdir() do path
        unshift!(LOAD_PATH, path)
        file_name = joinpath(path, string(module_name, ".jl"))
        sleep(2)
        touch(file_name)
        code = """module $(module_name)\nend\n"""
        write(file_name, code)
        reload(module_name)
        @test isa(root_module(Symbol(module_name)), Module)
        @test shift!(LOAD_PATH) == path
        rm(file_name)
    end
end

# Issue #19960
let
    test_workers = addprocs(1)
    push!(test_workers, myid())
    save_cwd = pwd()
    temp_path = mktempdir()
    try
        cd(temp_path)
        load_path = mktempdir(temp_path)
        load_cache_path = mktempdir(temp_path)

        ModuleA = :Issue19960A
        ModuleB = :Issue19960B

        write(joinpath(load_path, "$ModuleA.jl"),
            """
            __precompile__(true)
            module $ModuleA
                export f
                f() = myid()
            end
            """)

        write(joinpath(load_path, "$ModuleB.jl"),
            """
            __precompile__(true)
            module $ModuleB
                using $ModuleA
                export g
                g() = f()
            end
            """)

        @everywhere test_workers begin
            unshift!(LOAD_PATH, $load_path)
            unshift!(Base.LOAD_CACHE_PATH, $load_cache_path)
        end
        try
            @eval using $ModuleB
            uuid = Base.module_uuid(root_module(ModuleB))
            for wid in test_workers
                @test Base.Distributed.remotecall_eval(Main, wid, :( Base.module_uuid(Base.root_module($(QuoteNode(ModuleB)))) )) == uuid
                if wid != myid() # avoid world-age errors on the local proc
                    @test remotecall_fetch(g, wid) == wid
                end
            end
        finally
            @everywhere test_workers begin
                shift!(LOAD_PATH)
                shift!(Base.LOAD_CACHE_PATH)
            end
        end
    finally
        cd(save_cwd)
        rm(temp_path, recursive=true)
        pop!(test_workers) # remove myid
        rmprocs(test_workers)
    end
end

end # !withenv

max_stars_repo_name: goretkin/julia

max_issues_repo_name: goretkin/julia

max_forks_repo_name: goretkin/julia


text:
# This file is a part of ChangesOfVariables.jl, licensed under the MIT License (MIT).

using ChangesOfVariables
using Test

using LinearAlgebra

using ChangesOfVariables: test_with_logabsdet_jacobian

include("getjacobian.jl")

@testset "torv_and_back" begin
    for x in (rand(3), 0.5, Complex(0.2,0.7), (3,5,9), Ref(42), rand(3, 4, 5), Complex.(rand(3,5), rand(3,5)))
        V, to_x = torv_and_back(x)
        @test V isa AbstractVector{<:Real}
        @test V == torv_and_back(x)[1]
        @test x isa Ref ? to_x(V)[] == x[] : to_x(V) == x
    end
end

@testset "test_with_logabsdet_jacobian" begin
    rx = 0.5
    cx = Complex(0.2, -0.7)
    X = rand(3, 3)
    CX = Complex.(randn(3,3), randn(3,3))

    myisapprox(a, b; kwargs...) = isapprox(a, b; kwargs...)

    noninferrable_inv(x) = x!=rand(size(x)...) ? inv(x) : ""
    ChangesOfVariables.with_logabsdet_jacobian(::typeof(noninferrable_inv), x) = noninferrable_inv(x), with_logabsdet_jacobian(inv, x)[2]
    @test_throws ErrorException @inferred with_logabsdet_jacobian(noninferrable_inv, rand(2, 2))

    test_with_logabsdet_jacobian(inv, rx, ForwardDiff.derivative, atol = 10^-6)
    test_with_logabsdet_jacobian(inv, cx, getjacobian, atol = 10^-6)
    test_with_logabsdet_jacobian(inv, X, ForwardDiff.jacobian, atol = 10^-6)
    test_with_logabsdet_jacobian(inv, CX, getjacobian, atol = 10^-6)
    test_with_logabsdet_jacobian(inv, CX, getjacobian, atol = 10^-6)
    test_with_logabsdet_jacobian(inv, CX, getjacobian, compare = myisapprox, atol = 10^-6)
    test_with_logabsdet_jacobian(noninferrable_inv, CX, getjacobian, atol = 10^-6)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SeqGen_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SeqGen")
JLLWrappers.@generate_main_file("SeqGen", UUID("1820de90-959d-54bb-9e2e-c539b020f79b"))
end  # module SeqGen_jll