Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_DistriCorrection.f90 @ 1750

Last change on this file since 1750 was 1609, checked in by lfita, 7 years ago

Adding the new fortrans:

DistriCorrection?.f90: Program to test the modification of the distribution
module_definitions.f90: Module with the general definitions
module_scientific.f90: Module with the scientific calculations
module_DistriCorection.f90: Module which modifies a distribution of values

File size: 15.6 KB

Line
1	MODULE module_DistriCorrection
2	! Modules for modIFication of a Ditribution
3
4	USE module_definitions
5	USE module_generic
6	USE module_scientific
7
8	CONTAINS
9
10	!!!!!!! Subroutines
11	! randLowRemovePercent_HighEqui: Method to re-distribute from lower tail of thedistribution to the
12	! higher one following percentages preserving total amount.
13
14	SUBROUTINE randLowRemovePercent_HighEqui3D(d1, d2, d3, dist, drng, rangeVals, dpercen, PercentRange, &
15	KindRemove, Nquants, Hquant, RepvalS, fillValue, dbg, newdist, redistributed, rngamount, &
16	Habovequant, aboveind)
17	! Method to re-distribute from lower tail of thedistribution to the higher one following percentages
18	! preserving total amount. This method, removes values from the lower end of the distribution and
19	! equi-distribute the total amount on a percentile on the high end
20
21	IMPLICIT NONE
22
23	INTEGER, INTENT(in) :: d1, d2, d3, drng, dpercen, Nquants
24	REAL(r_k), INTENT(in) :: Hquant, fillValue
25	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: dist
26	REAL(r_k), DIMENSION(drng), INTENT(in) :: rangeVals
27	REAL(r_k), DIMENSION(dpercen), INTENT(in) :: PercentRange
28	CHARACTER(len=50), INTENT(in) :: KindRemove, RepvalS
29	LOGICAL, INTENT(in) :: dbg
30	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: newdist
31	REAL(r_k), DIMENSION(d1,d2,drng), INTENT(out) :: rngamount
32	INTEGER, DIMENSION(d1,d2,d3), INTENT(out) :: redistributed, aboveind
33	REAL(r_k), DIMENSION(d1,d2), INTENT(out) :: Habovequant
34
35	! Local
36	INTEGER :: i, j, iv, ir, irm, ind, ia
37	INTEGER :: ierr
38	INTEGER :: irg, erg
39	INTEGER :: iHquant
40	INTEGER :: Nvals, Nrngvals, Nrmvals, Nredis, Nrmprev
41	INTEGER :: NaboveHquants
42	INTEGER, DIMENSION(d3) :: distind
43	REAL(r_k), DIMENSION(d3) :: varvals
44	REAL(r_k), DIMENSION(:),ALLOCATABLE :: varvalues, quants, newquants
45	REAL(r_k) :: one, zero
46	REAL(r_k) :: ndist, xdist, meandist
47	REAL(r_k) :: Repval, leakvals
48	REAL(r_k) :: iamount, initial_amount, eamount
49	REAL(r_k) :: TOTamount, dTOTamount
50	INTEGER, DIMENSION(:), ALLOCATABLE :: rngind, rmind
51	CHARACTER(len=4) :: IS1, IS2
52	CHARACTER(len=50) :: spc, RS1
53	CHARACTER(len=50) :: fname
54	CHARACTER(len=1000) :: vectorS1, vectorS2, bspc
55	LOGICAL, DIMENSION(d3) :: search1, search2, search
56	CHARACTER(len=256) :: msg
57
58	!!!!!!! Variables
59	! d[1/2/3]: size of the values to treat
60	! dist: values to re-distribute
61	! drng: number of ranges provided
62	! rangeVals: ',' values from which establish the ranges of removals
63	! dpercen: number of percentages provided
64	! PercentRange: ',' percentages of removing to apply to each range (one value less than [rangeVals])
65	! KindRemove: kind of methodology to select which values to remove
66	! 'rand': randomly selection within the correspondent range
67	! Nquants: quantity of percentiles to compute
68	! Hquant: percentile of the distribution of values above which total removed amount will be equi-distributed
69	! RepvalS: replacing value to be assigned at the removed values
70	! 'min': initial minimal value of the distribution
71	! 'mean': initial mean value of the distribution
72	! '[value]': any given value otherwise
73	! fillValue: missing value
74	!
75	! Outcome:
76	! newdist: Values of the new distribution
77	! redistributed: List of positions changed for each range
78	! rngamount: Amount re-distributed from each range
79	! Habovequant: value from which the amount is distributed
80	! aboveind: Indices of the values which received the distributed amount
81
82	fname = 'randLowRemovePercent_HighEqui'
83
84	spc = ' '
85	bspc = Nstrings(spc,4)
86	one = 1.
87	zero = 0.
88
89	IF (dbg) PRINT *," Re-aranging distribution of ", d3," values following '"//TRIM(fname)//"' method"
90
91	! Checking for the right amount of ranges and percentages
92	IF (dpercen /= drng - 1) THEN
93	WRITE(IS1,'(I4)')drng
94	WRITE(IS2,'(I4)')dpercen
95	vectorS1 = vectorR_KS(dpercen, PercentRange)
96	vectorS2 = vectorR_KS(drng, rangeVals)
97
98	msg = 'provided number of percentages: ' // TRIM(vectorS1) // ' (' // TRIM(IS2) // &
99	') does not fit for: ' // TRIM(IS1) // 'ranges:' // TRIM(vectorS2)
100	CALL ErrMsg(msg, fname, -1)
101	END IF
102
103	IF (dbg) THEN
104	PRINT *,TRIM(bspc) // 'Ranges modIFication _______'
105	DO ir=1, drng-1
106	WRITE(RS1,'(F5.2)')PercentRange(ir)*100
107	PRINT *,TRIM(bspc) // ' ' // TRIM(RS1) // ' % of (', rangeVals(ir), ',', rangeVals(ir+1), ']'
108	END DO
109	END IF
110
111	IF (ALLOCATED(quants)) DEALLOCATE(quants)
112	ALLOCATE(quants(Nquants+1), STAT=ierr)
113	msg = "Problem allocating 'quants'"
114	CALL ErrMsg(msg, fname, ierr)
115
116	! Checking for the upper-quantile
117	quants = RangeR_K(Nquants+1, zeroRK, 100.*oneRK)
118
119	iHquant = Index1DArrayR_K(quants, Nquants+1, hquant100.oneRK)
120
121	IF (iHquant == -1) THEN
122	WRITE(RS1,'(F50.10)')hquant*100.
123	WRITE(IS1,'(I4)')Nquants
124	vectorS1 = vectorR_KS(Nquants+1, quants)
125	msg = 'provided high percentile ' // TRIM(RS1) // ' does not fit with ' // TRIM(IS1) // &
126	' quantiles: ' // TRIM(vectorS1)
127	CALL ErrMsg(msg, fname, iHquant)
128	END IF
129
130	IF (dbg) PRINT *,TRIM(bspc) // 'To be added above:', quants(iHquant), ' [', iHquant, '] percentile'
131
132	! getting the right replacing value
133	SELECT CASE (TRIM(repvalS))
134	CASE ('min')
135	repval = ndist
136	CASE ('mean')
137	repval = meandist
138	CASE DEFAULT
139	READ(RepvalS,'(F50.25)', IOSTAT=ierr)Repval
140	msg = "For not 'min' or 'mean', a value is required as RepValS" // CHAR(10) // " provided: '" //&
141	TRIM(RepvalS) // "'"
142	CALL ErrMsg(msg, fname, ierr)
143	END SELECT
144
145	! Future new distribution
146	newdist = dist
147
148	! Total amounts for each range
149	rngamount = zero
150
151	! Distribution indices
152	distind = rangeI(d3,1,d3)
153
154	! Initialization of indices
155	redistributed = -1
156	aboveind = -1
157
158	!! Looping on all dimensions
159	DO i=1, d1
160	DO j=1, d2
161
162	varvals = dist(i,j,:)
163	search = (varvals == fillValue)
164	Nvals = COUNT(.NOT.search)
165
166	doingnothing: IF (Nvals == 0) THEN
167	PRINT *,TRIM(infomsg)
168	PRINT *,' ' // TRIM(fname) // ': all missed values for point:', i,', ',j, '!!'
169	newdist(i,j,:) = dist(i,j,:)
170	CYCLE
171
172	ELSE
173
174	IF (ALLOCATED(varvalues)) DEALLOCATE(varvalues)
175	ALLOCATE(varvalues(Nvals))
176
177	varvalues = PACK(varvals, MASK=.NOT.search)
178
179	! Amount of previous indices removed
180	Nrmprev = 0
181
182	ndist = MINVAL(varvalues)
183	xdist = MAXVAL(varvalues)
184	meandist = SUM(varvalues)/(d3*one)
185
186	IF (xdist < rangeVals(drng)) THEN
187	PRINT *,TRIM(infomsg)
188	PRINT *,' ' // TRIM(fname) // ': all values for point:', i,', ',j, ' are below ' // &
189	'ranges upper-end:', rangeVals(drng), ' !!'
190	newdist(i,j,:) = dist(i,j,:)
191	CYCLE
192	END IF
193
194	CALL quantilesR_K(d3, varvalues, Nquants+1, quants)
195
196	iamount = SUM(varvalues)
197	initial_amount = iamount + 0.
198
199	IF (dbg) THEN
200	PRINT *,TRIM(bspc) // 'values _______'
201	PRINT *,TRIM(bspc) // ' initial total amount:', iamount
202	PRINT *,TRIM(bspc) // ' min:', ndist, 'max:', xdist, 'mean:', meandist
203	PRINT *,TRIM(bspc) // ' quantiles:'
204	CALL PrintQuantilesR_K(d3, varvalues, Nquants+1, quants, Nstrings(spc,6))
205	END IF
206
207	IF (dbg) THEN
208	PRINT *, TRIM(bspc) // 'Loooping ranges _______'
209	END IF
210	bspc = Nstrings(spc,6)
211
212	! Getting total amount to remove and remove it from the value
213	TOTamount = zero
214
215	ranges: DO ir=1, drng-1
216	irg = ir
217	erg = ir + 1
218
219	IF (dbg) THEN
220	WRITE(RS1, '(F6.2)')PercentRange(ir)*100
221	PRINT ,TRIM(bspc) // ' re-aranging ' // TRIM(RS1) // ' % of interval (', rangeVals(irg), ', ', &
222	rangeVals(erg), ']'
223	END IF
224
225	! Searching the indices
226	WHERE (varvalues > rangeVals(irg))
227	search1 = .TRUE.
228	ELSEWHERE
229	search1 = .FALSE.
230	END WHERE
231
232	WHERE (varvalues <= rangeVals(erg))
233	search2 = .TRUE.
234	ELSEWHERE
235	search2 = .FALSE.
236	END WHERE
237
238	WHERE (search1 .AND. search2)
239	search = .TRUE.
240	ELSEWHERE
241	search = .FALSE.
242	END WHERE
243
244	! Values within the range
245	Nrngvals = COUNT(search)
246	IF (ALLOCATED(rngind)) DEALLOCATE(rngind)
247	ALLOCATE(rngind(Nrngvals), STAT=ierr)
248	msg = "Problem allocating 'rngind'"
249	CALL ErrMsg(msg, fname, ierr)
250
251	irm = 1
252	DO iv=1, d3
253	IF (search(iv)) THEN
254	rngind(irm) = iv
255	irm = irm + 1
256	END IF
257	END DO
258
259	! Number of values to remove
260	Nrmvals = INT(Nrngvals*PercentRange(ir))
261
262	IF (Nrmvals < 2) THEN
263	PRINT *,TRIM(warnmsg)
264	PRINT *,' ' // TRIM(fname) // ': very few values to remove:', Nrmvals, ' !!'
265	PRINT *,' for '//TRiM(RS1)//' % of interval (', rangeVals(irg), ',', rangeVals(erg), ']'
266	END IF
267
268	! Indices to remove
269	IF (ALLOCATED(rmind)) DEALLOCATE(rmind)
270	ALLOCATE(rmind(Nrmvals), STAT=ierr)
271	msg = "Problem allocating 'rmind'"
272	CALL ErrMsg(msg, fname, ierr)
273
274	! Removing kind
275	Kremove: SELECT CASE (TRIM(KindRemove))
276
277	CASE ('rand')
278	CALL rand_sample(Nrngvals, Nrmvals, rmind)
279
280	CASE DEFAULT
281	msg = "Removing kind '" // TRIM(KindRemove) // "' not ready" // CHAR(10) // &
282	"Available ones: 'rand'"
283	CALL ErrMsg(msg, fname, -1)
284
285	END SELECT Kremove
286
287	leakvals = zero
288	DO irm=1, Nrmvals
289	ind = rngind(rmind(irm))
290	! keeping the indices of the redistributed data
291	redistributed(i,j,Nrmprev+irm) = ind
292	rngamount(i,j,ir) = rngamount(i,j,ir) + (varvalues(ind) - repval)
293	! Keeping track of what is lost or gain during the transformation (for repval /= 0.)
294	leakvals = leakvals + (repval - varvalues(ind))
295
296	newdist(i,j,ind) = repval
297	END DO
298
299	TOTamount = TOTamount + rngamount(i,j,ir)
300	IF (dbg) THEN
301	PRINT *,TRIM(bspc) // ' removing:', Nrmvals, 'values from:', Nrngvals
302	PRINT *,TRIM(bspc) // ' from indices:', rngind, 're-distributing:', rngind(rmind(1:Nrmvals))
303	PRINT *,TRIM(bspc) // ' total amount to remove:', TOTamount
304	END IF
305	Nrmprev = Nrmprev + Nrmvals
306
307	END DO ranges
308
309	NonullRM: IF (TOTamount /= zero) THEN
310	! Equi-distributing the TOTamount on values >= hquant
311	Habovequant(i,j) = quants(iHquant)
312
313	IF (Habovequant(i,j) < rangeVals(drng)) THEN
314	PRINT *,TRIM(infomsg)
315	PRINT *,' point:', i,', ',j,' has a High-quantile:', quants(iHquant), ' smaller' // &
316	' than upper limit of ranges:', rangeVals(drng), '!!'
317	!PRINT *,' Trying to re-distribute on all the values above:', rangeVals(drng)
318	Habovequant(i,j) = oneRK
319
320	! Searching the indices
321	!WHERE (varvalues >= rangeVals(drng))
322	WHERE (varvalues > oneRK)
323	search = .TRUE.
324	ELSEWHERE
325	search = .FALSE.
326	END WHERE
327	PRINT *,' Re-distribute ', TOTamount, ' on', COUNT(search),' values above:', oneRK
328
329	ELSE
330	! Searching the indices
331	WHERE (varvalues >= Habovequant(i,j))
332	search = .TRUE.
333	ELSEWHERE
334	search = .FALSE.
335	END WHERE
336	END IF
337
338	! Quantity to add at each above quantile
339	NaboveHquants = COUNT(search)
340	dTOTamount = TOTamount / NaboveHquants
341
342	IF (dbg) THEN
343	PRINT *,TRIM(bspc) // 'incrementing:', NaboveHquants, 'values above:', Habovequant(i,j), &
344	' by', dTOTamount
345	PRINT *,TRIM(bspc) // ' replacing value:', repval
346	END IF
347
348	irm = 1
349	DO iv=1, d3
350	IF (search(iv)) THEN
351	aboveind(i,j,irm) = iv
352	newdist(i,j,iv) = varvalues(iv) + dTOTamount
353	irm = irm + 1
354	END IF
355	END DO
356
357	! Final total amount of the new distribution
358	eamount = SUM(newdist(i,j,:))
359
360	IF (dbg) THEN
361	Nredis = COUNT(redistributed(i,j,:) /= -1)
362	PRINT *,TRIM(bspc), Nredis, 'Values removed from ranges ________'
363	DO ia=1, Nredis
364	ind = redistributed(i,j,ia)
365	PRINT *,TRIM(bspc) // ' ', ia, ' prev:', varvalues(ind), 'new:', newdist(i,j,ind)
366	END DO
367	PRINT *,TRIM(bspc), NaboveHquants, 'Values above ________'
368	DO ia=1, NaboveHquants
369	ind = aboveind(i,j,ia)
370	PRINT *,TRIM(bspc) // ' ', ia, 'prev:', varvalues(ind), 'new:', newdist(i,j,ind)
371	END DO
372	END IF
373
374	! Check
375	IF (ABS(iamount - eamount) > 1.e-5) THEN
376	PRINT *, warnmsg
377	PRINT *,' ' // TRIM(fname) // ': final quantity:', eamount, ' &', iamount, ' does not' // &
378	' fit: abs(iamount - eamount) > 1.e-5:', abs(iamount-eamount), ' TOTamount:', TOTamount
379	STOP
380	END IF
381
382	IF (dbg) THEN
383	varvalues = newdist(i,j,:)
384	ndist = MINVAL(varvalues)
385	xdist = MAXVAL(varvalues)
386	meandist = SUM(varvalues)/(d3*one)
387
388	IF (ALLOCATED(newquants)) DEALLOCATE(newquants)
389	ALLOCATE(newquants(Nquants+1), STAT=ierr)
390	msg = "Problem allocating 'newquants'"
391	CALL ErrMsg(msg, fname, ierr)
392
393	CALL quantilesR_K(d3, varvalues, Nquants+1, newquants)
394
395	PRINT *,TRIM(bspc) // ' new distribution values _______'
396	PRINT *,TRIM(bspc) // ' final total amount:', eamount
397	PRINT *,TRIM(bspc) // ' min:', ndist, 'max:', xdist, 'mean:', meandist
398	PRINT *,TRIM(bspc) // ' initial quantiles:'
399	CALL PrintQuantilesR_K(20, varvalues, Nquants+1, quants, Nstrings(spc,8))
400	PRINT *,TRIM(bspc) // ' new quantiles:'
401	CALL PrintQuantilesR_K(20, varvalues, Nquants+1, newquants, Nstrings(spc,8))
402	!STOP
403	END IF
404
405	ELSE
406
407	PRINT *,TRIM(warnmsg)
408	PRINT *,' ' // TRIM(fname) // ': Any value has been removed !!'
409	Habovequant(i,j) = zero
410
411	END IF NonullRM
412
413	END IF doingnothing
414	END DO
415	PRINT *,' ' // TRIM(fname) // ': done', i, ' of', d1
416	END DO
417
418	IF (ALLOCATED(varvalues)) DEALLOCATE(varvalues)
419	IF (ALLOCATED(rngind)) DEALLOCATE(rngind)
420	IF (ALLOCATED(rmind)) DEALLOCATE(rmind)
421	IF (ALLOCATED(quants)) DEALLOCATE(quants)
422	IF (ALLOCATED(newquants)) DEALLOCATE(newquants)
423
424	RETURN
425
426	END SUBROUTINE randLowRemovePercent_HighEqui3D
427
428	END MODULE module_DistriCorrection

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format